cbm-rbd-paper.tex

\documentclass[authoryear]{myelsartcl}

% ------------ packages -------------

\usepackage[utf8]{inputenc}
\usepackage[OT1]{fontenc}
\usepackage{graphicx}
\usepackage[english]{babel}

\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{bm}

\usepackage[usenames,dvipsnames]{xcolor}
\usepackage{booktabs}
\usepackage{tikz}
\usepackage{sidecap}

\usepackage{url}
\usepackage[bookmarks]{hyperref}

%\usetikzlibrary{shapes.misc,fit}
\usetikzlibrary{%
   arrows,%
   calc,%
   fit,%
   patterns,%
   plotmarks,%
   shapes.geometric,%
   shapes.misc,%
   shapes.symbols,%
   shapes.arrows,%
   shapes.callouts,%
   shapes.multipart,%
   shapes.gates.logic.US,%
   shapes.gates.logic.IEC,%
   er,%
   automata,%
   backgrounds,%
   chains,%
   topaths,%
   trees,%
   petri,%
   mindmap,%
   matrix,%
   calendar,%
   folding,%
   fadings,%
   through,%
   patterns,%
   positioning,%
   scopes,%
   decorations.fractals,%
   decorations.shapes,%
   decorations.text,%
   decorations.pathmorphing,%
   decorations.pathreplacing,%
   decorations.footprints,%
   decorations.markings,%
   shadows}

%\usepackage{hyperref}
%\usepackage[bookmarks]{hyperref}
%\usepackage[colorlinks=true,citecolor=red,linkcolor=black]{hyperref}

% ------------ custom defs -------------

\newcommand{\reals}{\mathbb{R}}
\newcommand{\posreals}{\reals_{>0}}
\newcommand{\posrealszero}{\reals_{\ge 0}}
\newcommand{\naturals}{\mathbb{N}}

\newcommand{\dd}{\,\mathrm{d}}

\newcommand{\mbf}[1]{\mathbf{#1}}
\newcommand{\bs}[1]{\boldsymbol{#1}}
\renewcommand{\vec}[1]{{\bm#1}}

\newcommand{\uz}{^{(0)}} % upper zero
\newcommand{\un}{^{(n)}} % upper n
\newcommand{\ui}{^{(i)}} % upper i

\newcommand{\ul}[1]{\underline{#1}}
\newcommand{\ol}[1]{\overline{#1}}

\newcommand{\Tsys}{T_\text{sys}}

\newcommand{\Rsys}{R_\text{sys}}
\newcommand{\lRsys}{\ul{R}_\text{sys}}
\newcommand{\uRsys}{\ol{R}_\text{sys}}

\newcommand{\fsys}{f_\text{sys}}
\newcommand{\Fsys}{F_\text{sys}}
\newcommand{\lFsys}{\ul{F}_\text{sys}}
\newcommand{\uFsys}{\ol{F}_\text{sys}}

\newcommand{\lgt}{\ul{g}}
\newcommand{\ugt}{\ol{g}}

\newcommand{\E}{\operatorname{E}}
\newcommand{\V}{\operatorname{Var}}
\newcommand{\wei}{\operatorname{Wei}} % Weibull Distribution
\newcommand{\ig}{\operatorname{IG}}   % Inverse Gamma Distribution

\newcommand{\El}{\ul{\operatorname{E}}}
\newcommand{\Eu}{\ol{\operatorname{E}}}

\def\yz{y\uz}
\def\yn{y\un}
%\def\yi{y\ui}
\newcommand{\yfun}[1]{y^{({#1})}}
\newcommand{\yfunl}[1]{\ul{y}^{({#1})}}
\newcommand{\yfunu}[1]{\ol{y}^{({#1})}}

\def\ykz{y\uz_k}
\def\ykn{y\un_k}

\def\yzl{\ul{y}\uz}
\def\yzu{\ol{y}\uz}
\def\ynl{\ul{y}\un}
\def\ynu{\ol{y}\un}
\def\yil{\ul{y}\ui}
\def\yiu{\ol{y}\ui}

\def\ykzl{\ul{y}\uz_k}
\def\ykzu{\ol{y}\uz_k}
\def\yknl{\ul{y}\un_k}
\def\yknu{\ol{y}\un_k}

\newcommand{\ykzfun}[1]{y\uz_{#1}}
\newcommand{\ykzlfun}[1]{\ul{y}\uz_{#1}}
\newcommand{\ykzufun}[1]{\ol{y}\uz_{#1}}


\def\nz{n\uz}
\def\nn{n\un}
%\def\ni{n\ui}
\newcommand{\nfun}[1]{n^{({#1})}}
\newcommand{\nfunl}[1]{\ul{n}^{({#1})}}
\newcommand{\nfunu}[1]{\ol{n}^{({#1})}}

\def\nkz{n\uz_k}
\def\nkn{n\un_k}
\newcommand{\nkzfun}[1]{n\uz_{#1}}
\newcommand{\nkzlfun}[1]{\ul{n}\uz_{#1}}
\newcommand{\nkzufun}[1]{\ol{n}\uz_{#1}}

\def\nzl{\ul{n}\uz}
\def\nzu{\ol{n}\uz}
\def\nnl{\ul{n}\un}
\def\nnu{\ol{n}\un}
\def\nil{\ul{n}\ui}
\def\niu{\ol{n}\ui}

\def\nkzl{\ul{n}\uz_k}
\def\nkzu{\ol{n}\uz_k}
\def\nknl{\ul{n}\un_k}
\def\nknu{\ol{n}\un_k}

\def\yknow{y_k^{(\tnow)}}
\def\nknow{n_k^{(\tnow)}}

\newcommand{\nk}{n_k}
\newcommand{\nkp}{n_k'}
\newcommand{\yk}{y_k}
\newcommand{\ykp}{y_k'}

\def\taut{\tau(\vec{t})}
\def\ttau{\tilde{\tau}}
\def\ttaut{\ttau(\vec{t})}
\def\tautk{\tau(\vec{t}_k)}

\def\MZ{\mathcal{M}\uz}
\def\MN{\mathcal{M}\un}

\def\MkZ{\mathcal{M}\uz_k}
\def\MkN{\mathcal{M}\un_k}

\def\PkZ{\Pi\uz_k}
\def\PkN{\Pi\un_k}
\newcommand{\PZi}[1]{\Pi\uz_{#1}}

\def\tnow{t_\text{now}}
\def\tpnow{t^+_\text{now}}

\newcommand{\Rsysnow}{R^{(t_\text{now})}_\text{sys}}
\newcommand{\Tsysnow}{T^{(t_\text{now})}_\text{sys}}
\newcommand{\tsysnow}{t^{(t_\text{now})}_\text{sys}}
\newcommand{\fsysnow}{f^{(t_\text{now})}_\text{sys}}
\def\eknow{e_k^{(\tnow)}}
\def\cknow{c_k^{(\tnow)}}
\def\vectknow{\vec{t}_k^{(\tnow)}}
\def\Phinow{\Phi^{(\tnow)}}
\newcommand{\gnow}{g^{(\tnow)}}
\newcommand{\tausnow}{\tau_*^{(\tnow)}}
\newcommand{\tprep}{\tau_{\text{prep}}}
\newcommand{\tthresh}{\tau_{\text{thresh}}}
\newcommand{\tstarnow}{t_*^{(\tnow)}}
\newcommand{\gstarnow}{g_*^{(\tnow)}}
\newcommand{\gtotalnow}{g_\text{total}^{(\tnow)}}
\newcommand{\esys}{e_\text{sys}}
\newcommand{\mrsys}{\bar{r}_\text{sys}}

\newcommand{\tausfun}[1]{\tau_*^{(#1)}}


% ------------ options -------------

\allowdisplaybreaks

\journal{RESS}

\begin{document}

% ------------ frontmatter -------------

\begin{frontmatter}
\title{Condition-Based Maintenance for Complex Systems\\ based on Current Component Status\\ and Bayesian Updating of Component Reliability}

\author[tue]{Gero Walter}
\ead{g.m.walter@tue.nl}
\author[tue]{Simme Douwe Flapper}
\ead{s.d.p.flapper@tue.nl}

\address[tue]{School of Industrial Engineering, Eindhoven University of Technology, Eindhoven, The Netherlands}


\begin{abstract}
We propose a new condition-based maintenance policy for complex systems,
based on the status (working, defective) of all components within a system,
as well as the reliability block diagram of the system.
By means of the survival signature,
a generalization of the system signature allowing for multiple component types,
we obtain a predictive distribution for the system survival time,
also known as residual life distribution,
based on which of the system's components currently function or not,
and the current age of the functioning components.

The time to failure of the components of the system
is modeled by a Weibull distribution with a fixed shape parameter.
The scale parameter is iteratively updated in a Bayesian fashion
using the current (censored and non-censored) component lifetimes.
Each component type has a separate Weibull model that may also include test data.

The cost-optimal moment of replacement for the system is obtained by minimizing
the expected cost rate per unit of time.
The unit cost rate is recalculated when components fail
or at the end of every (very short) fixed inter-evaluation interval,
leading to a dynamic maintenance policy,
since the ageing of components and possible failures will change the cost-optimal moment of replacement in the course of time.
Via numerical experiments, some insight into the performance of the policy is given.
\end{abstract}

\begin{keyword}
condition-based maintenance \sep system reliability \sep remaining useful life \sep survival signature \sep unit time cost rate\end{keyword}
\end{frontmatter}


% ------------ manuscript -------------

\section{Introduction}
\label{intro}

Both in practice and academics,
there is a growing interest in condition-based maintenance (CBM), see e.g.,
\citet{2017:oldekeizer}
The central idea behind CBM is to maintain systems or components at exactly the right time,
i.e., just before they fail, %but not too early,
in order to keep their reliability high and operating costs low,
using information about the actual condition of the systems or components.
In this context a trade-off is made between the risk of failure during operation
(which can lead to costly downtime: idle workforce, missed production, penalties, loss of reputation)
and the costs of premature maintenance (wasting potential component or system lifetime,
downtime cost, cost of executing unnecessary maintenance activities).

Among the reasons for the above interest are the importance of increasingly short and
reliable delivery times, and decreasing profit margins due to worldwide competition.
All kinds of technical improvements have made it possible % it is more / better / earlier possible
to estimate the condition of systems as a whole as well as their underlying components,
whereas the cost for these technologies is decreasing rapidly,
making the use of CBM policies even more feasible.

Most CBM policies are based on a directly observable, continuously measurable condition or degradation signal,
e.g., the amount of vibration in case of rotating equipment.
Alternatively, such a signal or health status is constructed
using indirect measurements to determine the remaining useful life or time to failure of a system,
see, e.g., \citet{2014:rul-review, 2011:rul-review-statistical}.

In this paper, we propose a new kind of CBM policy
for when no such \emph{usual} degradation signal for the system is available,
but where the status (working or not working) of the system's components can be monitored (quasi-) continuously.
%i.e., observed every $\delta$ time units, where $\delta$ is very small compared to the lifetimes of components. 
In this situation, one can use the system's reliability block diagram
and information with respect to the status of its components
to directly calculate the system residual life distribution (RLD),
and base the maintenance policy for the system on this distribution.
In fact, our CBM policy can also be seen as based on a multivariate degradation signal,
where each component sends a binary signal, and the reliability block diagram is used for sensor fusion.

To calculate the system RLD, we use the survival signature,
which allows to deal with different types of components each having an individual failure time model
(see Section~\ref{sec:sysrel}).
Here we assume that components of different types fail independently,
and that failure times of components of the same type are independent given the model parameters.

In our new policy, there are three types of triggers to review the system:
(1) at the end of each very short time interval $\delta$ 
(where $\delta$ is very short compared to the lifetimes of components),
(2) directly after the failure of a component where the system still functions,
and (3) directly after the failure of a component due to which the system as a whole no longer can function. 
We denote each review of the system in reaction to triggers (1) and (2) by `evaluation',
since the system RLD is re-evaluated at each of these time points,
and a decision is made on whether preventive maintenance should be carried out or not.
We use the term `evaluation' in place of the term `inspection',
because the re-evaluation of the system RLD does not require a physical inspection,
as we assume to learn of the exact failure times of components online via (quasi-) continuous monitoring.

In case of a trigger of type (1),
the currently used models for component failure times are updated
to account for the extra lifetime observed for the components.
This information is taken into account via right-censored observations
in a Bayesian parameter update step.
In case of a trigger of type (2),
where a component has failed but the system is still functioning,
in addition to the component models also the reliability block diagram of the system is updated.
Based on the above changes,
it is determined in both cases whether the system should be replaced now,
or that this decision should be postponed until the next planned evaluation moment,
accepting the risk of a system breakdown before.
To determine the latter,
we calculate the `optimal' next evaluation moment %moment for system replacement***
based on minimizing the expected maintenance-related cost per unit of time for the present operational cycle.
An operational cycle spans the time between two consecutive system maintenance actions.
We assume that in both preventive and corrective maintenance, the entire system is replaced or brought to an as good as new state,
and thus optimize the unit cost rate calculated over the time since the last system replacement.

When the `optimal' next evaluation moment is before the next planned evaluation moment, we replace the system now;
else we wait for the next planned evaluation moment.
The reason behind this decision is that $\delta$ is taken as %usually is 
%the shortest interval for which re-evaluation of the system RLD will show any difference***
%the shortest inter-observation period that can be realized from a hardware point of view.
the shortest inter-evaluation period that can be realized from a hardware point of view. 
%the time granularity according to which maintenance actions can be reasonably triggered***.
To illustrate, for components with lifetimes amounting to a few days,
$\delta$ could be in the order of minutes;
for components with lifetimes measuring in months,
$\delta$ could be one day.

Note that strictly speaking, this decision procedure may not be optimal.
However, because $\delta$ is very small when compared to the prior estimates of the lifetime of the components,
at the time that the next `optimal' evaluation moment is before the next planned evaluation moment,
we usually already have done much better than without updating as described in this paper.

In this context it is also important to realize that our method
doesn't necessarily result in the actual optimal replacement solution for another reason,
inherent to using updated information for decisions,
when these updates can indicate that the former estimates were too optimistic or too pessimistic.
The reason for this is that for some components,
a too pessimistic prior probability density function (pdf)
may lead to an increase in system reliability over time without any intervention.
Whether or not this will be the case, as well as until when this kind of increase can occur,
is uncertain at every moment in time.
As a sidenote,
this kind of uncertainty about the probabilistic behaviour itself is very often neglected in operations research literature.
Determining the next `optimal' evaluation moment can thus only indicate whether or not waiting with the replacement
may result in lower or higher cost per unit of time,
i.e., only indicates a local optimum,
but not necessarily a global one.

In case of a trigger of type (3),
where the failure of a component leads to the whole system failing,
the currently used component failure time distributions are updated
like in case of a trigger of type (1) and (2)
to account for the information gained on components
until the moment of system failure.
Hereafter, a new system is installed,
using the updated models for all components based on the insight obtained from the finished operational cycles so far.

At first sight, the above quasi-continous updating of the component failure models,
and the system reliability block diagram,
as well as redetermining the best moment to replace the system,
seems to result in a (very) nervous maintenance policy.
However, due to the very short calculation times,
only at certain moments in time people involved in the execution of maintenance tasks
need to receive detailed information,
whereas during the remaining time,
all can be handled automatically, as currently done within, e.g., many safety systems. 

The setup of the rest of the paper is as follows.
In Section~\ref{sec:literature}, we describe our contribution to the CBM literature.
Next, in Sections~\ref{sec:sysrel}, \ref{sec:adaptive-sysrel-weibull} and \ref{sec:policy},
the three main steps in our policy are discussed in detail.
Hereafter, in Sections~\ref{sec:examples} and \ref{sec:sim},
our policy is applied further to an example.
Finally, in Section~\ref{sec:outlook},
a short summary and our main conclusions are given,
followed by some suggestions for further research.


\section{Literature review}
\label{sec:literature}

Although there is a lot of literature on CBM policies for multi-component systems,
see, e.g., \citet{2017:oldekeizer},
we found only a few papers having some relation to our approach.

The paper that comes closest to this paper is the paper by \citet{2013:si-et-al}.
In that paper, the condition of a system is inspected,
where the moment of the next inspection is determined via an updated degradation path based on the present condition of the system.
The authors provide exact expressions for the RLD,
also known as RUL (remaining useful lifetime) distribution,
which is updated in an empirical Bayesian framework using conjugate priors.
The RLD is used to construct a replacement decision model using a cost rate.

While we also use a Bayesian approach with conjugate priors to provide exact expressions for the RLD,
one of the main differences between \citet{2013:si-et-al} and our paper
is that we use component status data and the system reliability block diagram instead of a continuous degradation signal
as basis for calculating the system RLD.
Another main difference is that our policy can be applied to arbitrary system configurations
(i.e., arbitrary reliability block diagrams with arbitrary components),
whereas the method given by \citet{2013:si-et-al} can only be applied to single components.
A third difference concerns how the cost rate is calculated.
\citet{2013:si-et-al} calculate this rate using the renewal-reward theory approximation.
Essential for using the renewal-reward theory for calculating the expected cost per unit of time
is that the relevant characteristics of the system considered and its environment are not changing,
i.e.\ that the failure behavior of the system,
in essence its reliability function, does not change over time,
and that the same policy with the same parameter values is used for a statistically long enough period of time.
However, due to the update method related to the CBM policy presented in this paper, %like in \citet{2013:si-et-al},
we cannot use the renewal-reward theory based calculation method,
and therefore use the exact formula for the average cost per unit of time.

Another paper that includes some of the aspects covered in our paper is the paper by \citet{2012:sun-et-al}.
The authors construct a health index for a system based on sensor measurements,
where the health status prediction is updated sequentially,
resulting in a RUL distribution like with our method.
However, in contrast to our paper, the authors do not link their RUL distribution to any maintenance decisions.

\citet{2011:kim-et-al} developed a periodic monitoring CBM policy
where a maintenance decision is triggered when a Bayesian control chart
(a sequentially updated health indicator) exceeds a control limit (threshold)
that is determined by minimizing the expected average cost per time unit.
We use the same cost criterion, but base the maintenance policy decision directly on our exact RLD.

To summarize the main contribution of our paper,
it introduces a CBM-like policy that does not require knowledge of the degradation path of the system,
but can take into account continuous changes in the system reliability block diagram,
as well as changes in the component failure time distributions.



\section{Reliability function for complex systems using the survival signature}
\label{sec:sysrel}

Let $\Tsys$ be the time to failure of the system,
and $\Rsys(t) := P(\Tsys > t)$ be the system reliability function.
Here we describe how $\Rsys(t)$, given the system reliability block diagram and arbitrary component models,
can be efficiently calculated using the survival signature.
%We will show in Section~\ref{sec:adaptive-sysrel-weibull} how, for a certain choice of component model,
%this method can be used to find the system reliability function $\Rsysnow(t)$ at current time $\tnow$ for a monitored system,
%giving us the residual life distribution on which we base our adaptive maintenance policy.

We can analyze systems with an arbitrary reliability block diagram,
consisting of components of $K$ different types,
where there are $N_k$ exchangeable components of type $k$ in the system,
and $N$, the total number of components in the system, is given by $\sum_{k=1}^K N_k = N$.

As a running example, we consider a simplified automotive braking system
consisting of four component types $M$, $H$, $C$ and $P$,
with the reliability block diagram as depicted in Figure~\ref{fig:brakesys-layout}.
The master brake cylinder ($M$) activates all four wheel brake cylinders ($C_1$ -- $C_4$),
which in turn actuate a braking pad assembly each ($P_1$ -- $P_4$).
The hand brake mechanism ($H$) directly actuates the brake pad assemblies $P_3$ and $P_4$,
and the vehicle brakes when at least one brake pad assembly is actuated.
Note that because of the `handbrake shortcut', this system cannot be described as a nesting of series and parallel subsystems.
%This system is observed until time $\tnow$,
%leading to censored observation of certain component lifetimes within the system.
Our example system is deliberately kept quite simple for didactical reasons.
The survival signature allows to analyze much larger systems,
see, e.g., the examples in \citet{2017:reed},
and our approach is able to deal with such systems.

\begin{figure}
\centering
\begin{tikzpicture}
[typeM/.style={rectangle,draw,fill=black!20,thick,inner sep=0pt,minimum size=8mm}, %font=\footnotesize},
 typeC/.style={rectangle,draw,fill=black!20,thick,inner sep=0pt,minimum size=8mm}, %font=\footnotesize},
 typeP/.style={rectangle,draw,fill=black!20,thick,inner sep=0pt,minimum size=8mm}, %font=\footnotesize},
 typeH/.style={rectangle,draw,fill=black!20,thick,inner sep=0pt,minimum size=8mm}, %font=\footnotesize},
 type1/.style={rectangle,draw,fill=black!20,very thick,inner sep=0pt,minimum size=8mm},
 type2/.style={rectangle,draw,fill=black!20,very thick,inner sep=0pt,minimum size=8mm},
 type3/.style={rectangle,draw,fill=black!20,very thick,inner sep=0pt,minimum size=8mm},
 cross/.style={cross out,draw=red,very thick,minimum width=9mm, minimum height=7mm},
 hv path/.style={thick, to path={-| (\tikztotarget)}},
 vh path/.style={thick, to path={|- (\tikztotarget)}}]
\begin{scope}[xscale=1.5, yscale=1.2]
\node[typeM] (M)    at ( 0  , 0  ) {$M$};
\node[typeC] (C1)   at ( 1  , 1.5) {$C_1$};
\node[typeC] (C2)   at ( 1  , 0.5) {$C_2$};
%\node[cross]        at ( 1  , 0.5) {};
\node[typeC] (C3)   at ( 1  ,-0.5) {$C_3$};
%\node[cross]        at ( 1  ,-0.5) {};
\node[typeC] (C4)   at ( 1  ,-1.5) {$C_4$};
\node[typeP] (P1)   at ( 2  , 1.5) {$P_1$};
\node[typeP] (P2)   at ( 2  , 0.5) {$P_2$};
%\node[cross]        at ( 2  , 0.5) {};
\node[typeP] (P3)   at ( 2  ,-0.5) {$P_3$};
%\node[cross]        at ( 2  ,-0.5) {};
\node[typeP] (P4)   at ( 2  ,-1.5) {$P_4$};
\node[typeH] (H)    at ( 0  ,-1  ) {$H$};
\coordinate (start)  at (-0.7, 0);
\coordinate (startC) at ( 0.5, 0);
\coordinate (startH) at (-0.4, 0);
\coordinate (Hhop1)  at ( 0.4,-1);
\coordinate (Hhop2)  at ( 0.6,-1);
\coordinate (endP)   at ( 2.5, 0);
\coordinate (end)    at ( 2.8, 0);
\path (start)     edge[hv path] (M.west)
      (M.east)    edge[hv path] (startC)
      (startC)    edge[vh path] (C1.west)
                  edge[vh path] (C2.west)
                  edge[vh path] (C3.west)
                  edge[vh path] (C4.west)
      (C1.east)   edge[hv path] (P1.west)
      (C2.east)   edge[hv path] (P2.west)
      (C3.east)   edge[hv path] (P3.west)
      (C4.east)   edge[hv path] (P4.west)
      (endP)      edge[vh path] (P1.east)
                  edge[vh path] (P2.east)
                  edge[vh path] (P3.east)
                  edge[vh path] (P4.east)
                  edge[hv path] (end)
      (startH)    edge[vh path] (H.west)
      (H.east)    edge[hv path] (Hhop1)
      (Hhop1)     edge[thick,out=90,in=90] (Hhop2)
      (Hhop2)     edge[hv path] (P3.south)
                  edge[hv path] (P4.north);
\end{scope}
\end{tikzpicture}
\caption{Reliability block diagram for a simplified automotive brake system
with four component types $M$, $H$, $C$ and $P$.}
%Note that this system layout cannot be expressed as a nesting of series and parallel layouts.
%The corresponding survival signature $\Phi(l_M,l_H,l_C,l_P)$ is given in Table~\ref{tab:brakesys-survsign}.}
\label{fig:brakesys-layout}
\end{figure}

In a system with $N$ components, the state of the system can be expressed by the state vector
$\vec{x} = (x_1,x_2,\ldots,x_N) \in \{0,1\}^N$,
with $x_i=1$ if the $i$th component functions and $x_i=0$ if not.
The structure function $\phi : \{0,1\}^N \rightarrow \{0,1\}$, defined for all possible $\vec{x}$,
takes the value 1 if the system functions and 0 if the system does not function for state vector $\vec{x}$ \citep{BP75}.
Most real-life systems are coherent,
which means that $\phi(\vec{x})$ is non-decreasing in any of the components of $\vec{x}$,
so system functioning cannot be improved by worse performance of one or more of its components.
Furthermore, one can usually assume that $\phi(0, \ldots, 0) = 0$ and $\phi(1, \ldots, 1) = 1$.

The survival signature \citep{2012:survsign} is a summary of the structure function
for systems with $K$ groups of exchangeable components.
Denoted by $\Phi(l_1,\ldots,l_K)$, with $l_k=0,1,\ldots,N_k$ for $k=1,\ldots,K$,
it is defined as the probability for the event that the system functions
given that precisely $l_k$ of its $N_k$ components of type $k$ function, for each $k\in \{1,\ldots,K\}$.
Essentially, this creates a $K$-dimensional partition for the event $\Tsys > t$,
such that $\Rsys(t) = P(\Tsys > t)$ can be calculated using the law of total probability,
\begin{align}
P(\Tsys > t) &= \sum_{l_1=0}^{N_1} \cdots \sum_{l_K=0}^{N_K} P(\Tsys > t \mid C^1_t = l_1,\ldots, C^K_t = l_K)
                                                                                  P\Big( \bigcap_{k=1}^K \{ C^k_t = l_k\} \Big) \nonumber\\
             &= \sum_{l_1=0}^{N_1} \cdots \sum_{l_K=0}^{N_K} \Phi(l_1,\ldots,l_K) P\Big( \bigcap_{k=1}^K \{ C^k_t = l_k\} \Big) \nonumber\\
             &= \sum_{l_1=0}^{N_1} \cdots \sum_{l_K=0}^{N_K} \Phi(l_1,\ldots,l_K) \prod_{k=1}^K P(C^k_t = l_k)\,,
\label{eq:sysrel-survsign}
\end{align}
where $C^k_t$ is the number of type $k$ components functioning at time $t$,
and $P(C^k_t = l_k)$ is the (predictive) probability that exactly $l_k$ components of type $k$ function at time $t$.
The last equality holds as we assume that components of different types fail independently.
Note that for coherent systems, the survival signature $\Phi(l_1,\ldots,l_K)$ is non-decreasing in each $l_k$.

Continuing our example,
the survival signature for the system in Figure~\ref{fig:brakesys-layout} is given in Table~\ref{tab:brakesys-survsign},
omitting the entries for which $\Phi(l_M, l_H, l_C, l_P) = 0$ or $\Phi(l_M, l_H, l_C, l_P) = 1$,
since the full table would contain $\prod_{k=1}^K (N_k + 1) = 2 \times 2 \times 5 \times 5 = 100$ rows.
%
To illustrate how the table entries are obtained,
consider the case $l_M = 0, l_H = 1, l_C = 0, l_P = 1$,
i.e., $M$ has failed, $H$ functions, $C_1$ -- $C_4$ have failed, and one of $P_1$ -- $P_4$ functions.
The system functions only if the single functioning type $P$ component is either $P_3$ or $P_4$.
As the type $P$ components are exchangeable, each of $P_1$ -- $P_4$ has the same probability of being the functioning component,
and so the probability of $P_3$ or $P_4$ being the functioning component is $2 / 4$, thus $\Phi(0, 1, 0, 1) = 0.5$.
%
The omitted rows from Table~\ref{tab:brakesys-survsign} can be determined as follows:
\begin{itemize}
\item Take a row, decrease $l_k$ for one $k \in \{M, H, C, P\}$ by 1.
If the resulting row is not in Table~\ref{tab:brakesys-survsign}, then the corresponding value of $\Phi$ is $0$.
\item Take a row, increase $l_k$ for one $k \in \{M, H, C, P\}$ by 1.
If the resulting row is not in Table~\ref{tab:brakesys-survsign}, then the corresponding value of $\Phi$ is $1$.
\end{itemize} 
The \textsf{R} package \texttt{ReliabilityTheory} \citep{2016:aslett-RT}
provides a convenient function to calculate the survival signature table based on any given graph.


\begin{table}
\centering
\begin{tabular}{cccclcccccl}
  \toprule
$l_M$ & $l_H$ & $l_C$ & $l_P$ & $\Phi$ & \quad & $l_M$ & $l_H$ & $l_C$ & $l_P$ & $\Phi$\\ 
  \midrule
1 & 0 & 1 & 1 & 0.25 & & 1 & 0 & 2 & 1 & 0.50 \\ 
1 & 0 & 1 & 2 & 0.50 & & 1 & 0 & 2 & 2 & 0.83 \\ 
1 & 0 & 1 & 3 & 0.75 & & 1 & 0 & 3 & 1 & 0.75 \\ 
0 & 1 & 0 & 1 & 0.50 & & 1 & 1 & 0 & 1 & 0.50 \\ 
0 & 1 & 0 & 2 & 0.83 & & 1 & 1 & 0 & 2 & 0.83 \\ 
0 & 1 & 1 & 1 & 0.50 & & 1 & 1 & 1 & 1 & 0.62 \\ 
0 & 1 & 1 & 2 & 0.83 & & 1 & 1 & 1 & 2 & 0.92 \\ 
0 & 1 & 2 & 1 & 0.50 & & 1 & 1 & 2 & 1 & 0.75 \\ 
0 & 1 & 2 & 2 & 0.83 & & 1 & 1 & 2 & 2 & 0.97 \\ 
0 & 1 & 3 & 1 & 0.50 & & 1 & 1 & 3 & 1 & 0.88 \\ 
0 & 1 & 3 & 2 & 0.83 \\
0 & 1 & 4 & 1 & 0.50 \\
0 & 1 & 4 & 2 & 0.83 \\
   \bottomrule
\end{tabular}
\caption{Survival signature $\Phi(l_M, l_H, l_C, l_P)$
for the simplified automotive brake system depicted in Figure~\ref{fig:brakesys-layout},
omitting the rows for which $\Phi(l_M, l_H, l_C, l_P) = 0$ or $\Phi(l_M, l_H, l_C, l_P) = 1$.}
\label{tab:brakesys-survsign}
\end{table}


\section{Adaptive system residual life distribution based on Weibull component models}
\label{sec:adaptive-sysrel-weibull}

In this section, we describe how the system reliability function $\Rsysnow(t)$ at current time $\tnow$ 
can be calculated for a specific component model.
%Section~\ref{sec:policy} will then show how an adaptive maintenance policy
%can be derived for a monitored system based on such a current system residual life distribution. 

We consider the well-known Weibull model for the component lifetimes,
which is used in a wide variety of reliability studies. 
To keep things simple, we assume that the shape parameter of the Weibull distribution is known,
and that only the scale parameter needs to be estimated.
The model could be extended to learn also the shape parameter in a later step,
using, e.g., the discretized approach by \cite{1969:soland}.

By using the Bayesian approach described in Section~\ref{sec:weibull},
our component model allows to include both expert assessments and test data (if available),
and furthermore can account for component lifetime information from previous operational cycles
and the current operational cycle up to $\tnow$.
%***component model: weibull with fixed shape,
%inverse gamma prior on scale parameter,
%fix priors via expected lifetime and prior strength,
%test data inclusion
%with noninformatiove right-censoring
%
In Section~\ref{sec:postpred} we derive the posterior predictive distribution
for the Weibull component model as needed for the current system reliability calculation.
Finally, in Section~\ref{sec:sysreltnow}, we adapt the method from Section~\ref{sec:sysrel} to the dynamic setting,
resulting in a formula for $\Rsysnow(t)$.
%***output is $\Rsysnow(t)$, the current (at time $\tnow$) system reliability function
%(residual life distribution, RLD)
%taking into account the current system state,
%including the current ages of system components, 
%and the lifetime histories of all component types,
%including test data (if available) and expert assessments.

The model description below follows closely \cite[\S 2, \S 4.1 -- \S 4.3]{2016:walter-coolen},
who presented the same residual life distribution model.
However, \cite{2016:walter-coolen} do not consider maintenance policies,
but focus instead on a generalization of the RLD model by using imprecise probability models for the component lifetimes.
%These sets of priors can be seen as parametric P-boxes,
%a certain kind of imprecise probability model \citep[see, e.g.,][]{itip}.
%Such models allow for vague and partial prior specifications,
%and provide sensitivity to prior-data conflict \cite[\S 2.2.3.3]{diss}.


\subsection{Bayesian update of the Weibull component models}
\label{sec:weibull}

Here we describe the Weibull model for the component lifetimes,
together with the Bayesian update procedure which allows to include
expert knowledge, component tests, and information from the monitored components in the system.

For each type $k$ component, we assume that the lifetime $T_{k,i}$ ($i=1,\ldots,N_k$, $k = 1, \ldots, K$)
is Weibull distributed with scale parameter $\lambda_k$ and (fixed) shape parameter $\beta_k > 0$,
in short $T_{k,i} \mid \lambda_k \sim \wei(\beta_k,\lambda_k)$,
with pdf (probability density function) and cdf (cumulative distribution function)
\begin{align}
\label{eq:weibulldens}
f(t_{k,i} \mid \lambda_k) &= \frac{\beta_k}{\lambda_k} (t_{k,i})^{\beta_k-1} e^{-\frac{(t_{k,i})^{\beta_k}}{\lambda_k}}\,, \\
\label{eq:weibullcdf}
F(t_{k,i} \mid \lambda_k) &= 1 - e^{-\frac{(t_{k,i})^{\beta_k}}{\lambda_k}} = P(T_{k,i} \leq t_{k,i} \mid \lambda_k)\,,
\end{align}
where $\lambda_k > 0$ and $t > 0$.

The shape parameter $\beta_k$ determines whether the hazard rate is increasing ($\beta_k > 1$)
or decreasing ($\beta_k < 1$) over time.
For $\beta_k=1$, one obtains the Exponential distribution with constant hazard rate as a special case.
The scale parameter $\lambda_k$ can be interpreted through the relation
\begin{align}
\E[T_{k,i} \mid \lambda_k] &= \lambda_k^{1/\beta_k}\, \Gamma(1 + 1/\beta_k)\,.
\label{eq:lambdainterpret}
\end{align}
We will use this equation to convert expected lifetimes to $\lambda_k$ and vice versa.

With a Bayesian approach, one can express expert knowledge about the reliability of the components
by assigning a so-called prior distribution,
a distribution over the unknown parameter $\lambda_k$.
This prior distribution $f(\lambda_k)$ is then updated 
to the so-called posterior distribution $f(\lambda_k \mid \vec{t})$,
the distribution over $\lambda_k$ given the data $\vec{t}$,
via Bayes' rule
\begin{align*}
f(\lambda_k \mid \vec{t}) &\propto f(\vec{t}\mid\lambda_k) f(\lambda_k)\,.
\end{align*}
The posterior $f(\lambda_k \mid \vec{t})$ subsumes the information
from both expert knowledge and data,
and forms the basis for all inferences, like, e.g., predictions.

For the prior over $\lambda_k$,
a convenient choice is to use the inverse Gamma distribution,
which is commonly parametrized in terms of the parameters $a_k > 0$ and $b_k > 0$,
\begin{align}
f(\lambda_k\mid a_k,b_k) &= \frac{(b_k)^{a_k}}{\Gamma(a_k)} \lambda_k^{-a_k -1} e^{-\frac{b_k}{\lambda_k}}\,,
\label{eq:ig-def}
\end{align}
in short, $\lambda_k \mid a_k, b_k \sim \ig(a_k,b_k)$.
%Here, we have added the prior parameters $a_k$ and $b_k$ in notation
%to indicate that the prior over $\lambda$ depends on their values.

The inverse Gamma distribution is convenient because it is a conjugate prior,
i.e., the posterior obtained by Bayes' rule is again an inverse Gamma distribution and thus easily tractable
because only the prior parameters need to be updated to obtain the posterior parameters,
so no numerical integation or simulation techniques are necessary.
Furthermore, this conjugacy holds also when right-censored observations are used for updating,
as indicated below.

Instead of the usual parametrization in terms of $a_k$ and $b_k$,
we use the parameters $\nk > 1$ and $\yk > 0$
%(the upper index ${}\uz$ is used to indicate that these are prior parameters)
because they are more easy to interpret.
They are defined as
\begin{align}
\nk &= a_k - 1 & &\text{ and}
&
\yk &= b_k / \nk,
\label{eq:abtony}
\end{align}
where $\yk$ can be interpreted as the prior guess for the scale parameter $\lambda_k$,
as $\E[\lambda_k\mid\nk,\yk] = \yk$.
Using \eqref{eq:lambdainterpret},
we can thus translate an expert's statement of the expected component lifetime into a corresponding value for $\yk$.
$\nk$ can be seen as a prior strength or pseudocount,
this will become clear in the discussion of the update step below.

The parametrization in terms of $\nk$ and $\yk$ also clarifies the nature of the combination
of prior information and data through Bayes' rule.
In the conjugate setting,
applying Bayes' rule simply means that the prior parameters $\nk$ and $\yk$
are updated to posterior parameters, which we denote by $\nkp$ and $\ykp$, respectively.
Assume we observe $N_k = e_k + c_k$ component lifetimes,
where $e_k$ is the number of actual failure events,
and $c_k$ is the number of right-censored observations.
We denote the failure times by $t_{k,1}, \ldots, t_{k,e_k}$,
and the censoring times by $t^+_{k,1}, \ldots, t^+_{k,c_k}$,
and collect them in an observation vector $\vec{t}_k = (t_{k,1}, \ldots, t_{k,e_k}, t^+_{k,1}, \ldots, t^+_{k,c_k})$.
Then, the updated, posterior parameters are
\begin{align}
\nkp &= \nk + e_k\,, 
&
%\ykp &=  \frac{\nk \yk + \tautk}{\nk + e_k}\,,
\ykp &=  \frac{\nk}{\nk + e_k}\, \yk + \frac{e_k}{\nk + e_k}\, \frac{\tautk}{e_k}\,,
\label{eq:ig-update}
\end{align}
where $\tautk = \sum_{i=1}^{e_k} (t_{k,i})^\beta + \sum_{i=1}^{c_k} (t^+_{k,i})^\beta$.
%The upper index ${}\un$ indicates that these are posterior parameters resulting from an update with $n_k$ observations,
%where we leave out the index $k$ for $n_k$ in the upper index to increase legibility.

From the simple update rule \eqref{eq:ig-update}, we see that
$\ykp$ is a weighted average of the prior parameter $\yk$ and the maximum likelihood estimator (ML) $\tautk/e_k$,
with weights proportional to $\nk$ and $e_k$, respectively.
$\nk$ can thus be interpreted as a prior strength or pseudocount,
indicating how much our prior guess should weigh against the $e_k$ observed failure events.
Furthermore, $\V[\lambda\mid\nk,\yk] = (\yk)^2 / (1 - 1/\nk)$;
for fixed $\yk$, a higher $\nk$ indicates 
that more probability mass is concentrated around $\yk$.

Using \eqref{eq:abtony} and \eqref{eq:ig-update}, the posterior distribution over $\lambda_k$ is given by
\begin{align}
\lambda_k \mid \nk, \yk, \vec{t}_k \sim \ig(\nk + e_k + 1, \nk \yk + \tautk)\,.
\label{eq:ig-update-alpha}
\end{align}
%where we have added the prior parameters $\nkz$ and $\ykz$ in notation
%to emphasize that the posterior is a synthesis of prior information and data.
As this posterior can be defined in terms of
the updated, posterior parameters $\nkp$ and $\ykp$,
we may also write
\begin{align*}
f(\lambda_k \mid \nk, \yk, \vec{t}_k) &= f(\lambda_k \mid \nkp, \ykp)\,.
\end{align*}
The iterative nature of Bayesian inference means that we can take the updated,
posterior parameters as new prior parameters and update them again using new data.
We will use this property to use the extra information about components
that accumulates during operation of the system:
At any time $\tnow > 0$,
all non-failed components contribute a right censored observation $\tpnow$;
any failed components instead contribute a fully observed, non-censored lifetime,
and both types of observations can be used in the update step \eqref{eq:ig-update}.
%Before we write this out formally in Section~\ref{sec:postpred},
%we will first present how system reliability functions can be derived based on component models.

%To keep notation simple,
%we will below denote the parameters of the inverse Gamma distribution by $\nkz$ and $\ykz$,
%regardless of them expressing expert information alone,
%or stemming from the combination of expert information and test data.


\subsection{Component posterior predictive distributions determined at $\tnow$}
\label{sec:postpred}

To calculate the current system reliability function $\Rsysnow(t)$,
we need, for each $k=1,\ldots, K$, the probabilities $P(C^k_t = l_k)$
for the number of type $k$ components that function at times $t > \tnow$,
taking into account all information available at $\tnow$.
In the Bayesian framework, these probabilities are given
by a so-called posterior predictive distribution that can be derived from the posterior over $\lambda$.

Denote by $\nkz$ and $\ykz$ the parameters reflecting the knowledge base at system start-up time $t=0$.
These could have been obtained by updating prior parameters $\nk$ and $\yk$ (reflecting expert knowledge)
to $\nkp$ and $\ykp$ using test data,
or could be taken directly equal to $\nk$ and $\yk$ if no test data is available.
In both cases, the distribution $f(\lambda_k \mid \nkz, \ykz)$
thus accounts for all knowledge on component type $k$ that is available at system start-up.

Following our comments at the end of Section~\ref{sec:weibull},
the component models can be further updated
using information gained from the current operational cycle as of $\tnow$: %run of the system until $\tnow$:
failure times of failed components,
and the right-censoring time $\tpnow$ for each non-failed component,
can be used to update $\nkz$ and $\ykz$ according to \eqref{eq:ig-update}.
We denote the resulting parameter values by $\nknow$ and $\yknow$.
%(It is of course possible to update $\nkz$ and $\ykz$ directly to $\nknow$ and $\yknow$ if no test data is available.)
%Furthermore, the fact that the components in the system that still function have reached the age $\tnow$
%must be used when calculating the posterior predictive distribution.

In analogue to the notation used in Section~\ref{sec:weibull},
let $N_k = \eknow + \cknow$ be the number of type $k$ components in the monitored system,
where $\eknow$ denotes the number of type $k$ components that have failed by $\tnow$ during the current operational cycle,
and $\cknow$ denotes the number of type $k$ components that have not failed by $\tnow$ during the current operational cycle.
We can thus collect these observations, as of $\tnow$, in a vector
$\vectknow = (t_{k,1}, \ldots, t_{k,\eknow}, \tpnow, \ldots, \tpnow)$,
containing $\cknow$ right-censored observations $\tpnow$.
%
The resulting posterior predictive distribution for any $t > \tnow$ is obtained as
\begin{align}
\lefteqn{%
P(C^k_t = l_k\mid\nkz,\ykz, \vectknow) }\hspace*{5.75ex} \nonumber\\  %
 &= { \cknow \choose l_k} \int \big[P(T^k >    t \mid T^k > \tnow, \lambda_k)\big]^{l_k} \times \nonumber\\ & \hspace*{3.5ex}
                               \big[P(T^k \leq t \mid T^k > \tnow, \lambda_k)\big]^{\cknow - l_k}
    f(\lambda_k\mid\nkz,\ykz,\vectknow) \dd \lambda_k\,,
\label{eq:postpredtnow}
\end{align}
where $T^k$ is the Weibull distributed lifetime of a component of type $k$.
Note that through the condition $T^k > \tnow$, we also take into account
that the components in the system have the age $\tnow$.
Now, by the Weibull assumption \eqref{eq:weibullcdf}, one has
\begin{align}
P(T^k \leq t \mid T^k > \tnow, \lambda_k)
 &= \frac{P(\tnow < T^k \leq t \mid\lambda_k)}{P(T^k > \tnow \mid \lambda_k)} \nonumber\\
 &= \frac{F(t\mid\lambda_k) - F(\tnow\mid\lambda_k)}{1-F(\tnow\mid\lambda_k)} 
% = \frac{e^{-\frac{(\tnow)^{\beta_k}}{\lambda_k}} - e^{-\frac{t^{\beta_k}}{\lambda_k}}}{e^{-\frac{(\tnow)^{\beta_k}}{\lambda_k}}}
  = 1 - e^{-\frac{t^{\beta_k} - (\tnow)^{\beta_k}}{\lambda_k}}\,.
\label{eq:weibullcondprob}
\end{align}
Consider the posterior \eqref{eq:ig-update-alpha} obtained using observations $\vectknow$,
written in terms of the updated parameters $\nknow$ and $\yknow$.
Substituting this posterior and \eqref{eq:weibullcondprob} into \eqref{eq:postpredtnow}, we get
\begin{align}
\lefteqn{P(C^k_t = l_k\mid\nkz,\ykz, \vectknow)}\hspace*{5.75ex} \nonumber\\
 &= { \cknow \choose l_k} \int \Big[    e^{-\frac{t^{\beta_k} - (\tnow)^{\beta_k}}{\lambda_k}}\Big]^{l_k}
                               \Big[1 - e^{-\frac{t^{\beta_k} - (\tnow)^{\beta_k}}{\lambda_k}}\Big]^{\cknow - l_k}
    \times \nonumber\\ & \hspace*{14.5ex}
    \frac{\big(\nknow\yknow\big)^{\nknow + 1}}{\Gamma(\nknow + 1)}
    \lambda_k^{-(\nknow + 1) - 1} e^{-\frac{\nknow\yknow}{\lambda_k}} \dd \lambda_k \nonumber\\
 &= { \cknow \choose l_k} \sum_{j=0}^{\cknow - l_k} (-1)^j { \cknow - l_k \choose j}
    \frac{\big(\nknow\yknow\big)^{\nknow + 1}}{\Gamma(\nknow + 1)} 
    \times \nonumber\\ & \hspace*{1ex}
    \int \lambda_k^{-(\nknow + 1) - 1}
    \exp\Big\{-\frac{(l_k + j) (t^{\beta_k} - (\tnow)^{\beta_k}) + \nknow\yknow}{\lambda_k}\Big\} \dd \lambda_k\,.
\end{align}
The terms remaining under the integral form the core of an inverse Gamma distribution \eqref{eq:ig-def}
with parameters $\nknow + 1$ and $\nknow\yknow + (l_k + j) (t^{\beta_k} - (\tnow)^{\beta_k}))$,
allowing to solve the integral using the corresponding normalization constant.
We thus have, for $l_k \in \{0,1,\ldots,\cknow\}$,
\begin{align}
\label{eq:postpred-priorparams}
\lefteqn{P(C^k_t = l_k\mid\nkz,\ykz, \vectknow)}\hspace*{5ex} \nonumber\\
 &= { \cknow \choose l_k} \sum_{j=0}^{\cknow - l_k} (-1)^j { \cknow - l_k \choose j} \times \nonumber\\ & \hspace*{10ex}
    \left(\frac{\nknow\yknow}{\nknow\yknow + (l_k + j) \big(t^{\beta_k} - (\tnow)^{\beta_k}\big)}\right)^{\nknow + 1} \,, \\
 & \quad\text{where } \nknow\yknow = \nkz\ykz + \sum_{i=1}^{\eknow} (t_{k,i})^{\beta_k} + \cknow (\tnow)^{\beta_k} \nonumber\\
 & \quad\text{and }  \nknow = \nkz + \eknow \,. \nonumber
% &= { \cknow \choose l_k} \sum_{j=0}^{\cknow - l_k} (-1)^j { \cknow - l_k \choose j} \times \nonumber\\ & \hspace*{1ex}
%    \left(\frac{\nkz\ykz + \sum_{i=1}^{\eknow} (t_{k,i})^{\beta_k} + \cknow (\tnow)^{\beta_k} }%
%               {\nkz\ykz + \sum_{i=1}^{\eknow} (t_{k,i})^{\beta_k} + \cknow (\tnow)^{\beta_k} +%
%                (l_k + j) \big(t^{\beta_k} - (\tnow)^{\beta_k}\big) }\right)^{\nkz + \eknow + 1}.
% &= \sum_{j=0}^{\cknow - l_k} (-1)^j \frac{\cknow !}{l_k! j! (\cknow - l_k - j)!}   
%    \left(\frac{\nknow\yknow}{\nknow\yknow + (l_k + j) \big(t^{\beta_k} - (\tnow)^{\beta_k}\big)}\right)^{\nknow + 1} \nonumber\\
% &= \sum_{j=0}^{\cknow - l_k} (-1)^j \frac{\cknow !}{l_k! j! (\cknow - l_k - j)!} \times \nonumber\\ & \hspace*{10ex}  
%    \left(\frac{\nkz\ykz + \sum_{i=1}^{\eknow} (t_{k,i})^{\beta_k} + \cknow       (\tnow)^{\beta_k} }%
%               {\nkz\ykz + \sum_{i=1}^{\eknow} (t_{k,i})^{\beta_k} + (\cknow - l_k - j) (\tnow)^{\beta_k} + (l_k + j) t^{\beta_k} }\right)^{%
%    \nkz + \eknow + 1}.
\end{align}
%where $\nknow\yknow = \nkz\ykz + \sum_{i=1}^{\eknow} (t_{k,i})^{\beta_k} + \cknow (\tnow)^{\beta_k}$
%and $\nknow = \nkz + \eknow + 1$.
%These posterior predictive probabilities can also be expressed as a cumulative probability mass function (cmf)
%\begin{align}
%F(l_k \mid \nkz,\ykz,\vectknow) = P(C^k_t \leq l_k \mid \nkz,\ykz,\vectknow) 
% = \sum_{j=0}^{l_k} P(C^k_t = j \mid \nkz,\ykz,\vectknow)\,.
%\end{align}


\subsection{Conditional system reliability function determined at $\tnow$}
\label{sec:sysreltnow}

Now that we can calulate the component posterior predictive distributions at $\tnow$,
we can use the method from Section~\ref{sec:sysrel}
to determine the dynamic system residual life distribution (RLD) $\Rsysnow(t)$,
giving us the probability that the system functions at times $t > \tnow$,
taking into account all information available at $\tnow$:
\begin{align}
\Rsysnow(t) &= \sum_{l_1=0}^{c_1^{(\tnow)}} \cdots \sum_{l_K=0}^{c_K^{(\tnow)}} \Phinow(l_1,\ldots,l_K)
               \prod_{k=1}^K P(C^k_t = l_k\mid\nkz,\ykz, \vectknow)\,.
\label{eq:sysrel-tnow}
\end{align}
Since this is the distribution of system reliability conditional on the system surviving until $\tnow$,
we have that $\Rsysnow(\tnow) = 1$ for any $\tnow$ before system failure.
Note that if one or several components have failed by $\tnow$,
the system reliability block diagram changes, and with it the survival signature $\Phi(l_1,\ldots,l_K)$.
We denote the current survival signature by $\Phinow(l_1,\ldots,l_K)$.

For each prediction time $t$,
\eqref{eq:sysrel-tnow} is a sum over $\prod_{k=1}^K (\cknow + 1)$ terms.
However, some of these terms correspond to $\Phinow(l_1,\ldots,l_K) = 0$,
which can thus be disregarded.
For each of the remaining terms,
we must calculate the product $\prod_{k=1}^K P(C^k_t = l_k\mid\nkz,\ykz, \vectknow)$.
For all products, the constituting factors can be taken from the same table
enumerating $P(C^k_t = l_k\mid\nkz,\ykz, \vectknow)$ for all $l_k \in \{0, 1, \cknow\}$, $k=1\ldots,K$,
so \eqref{eq:postpred-priorparams} needs to be evaluated 
only $\sum_{k=1}^K (\cknow + 1)$ times.

We will evaluate \eqref{eq:sysrel-tnow} on a dense grid of prediction times $t > \tnow$,
thus discretely approximating the RLD.
As the evaluation of \eqref{eq:sysrel-tnow} for each $t$ does not involve any complex numeric calculations or Monte Carlo sampling,
the grid of prediction values can be very fine.
%
%Next, we will describe an adaptive maintenance policy based on %\eqref{eq:sysrel-tnow}
%a current residual life distribution,
%first in general, then for our implementation,
%where a fine discrete approximation of the RLD is obtained.


\section{Dynamic and adaptive maintenance policy based on the current system residual life distribution}
\label{sec:policy}

%***start with $\tausnow$???
Now we can derive the adaptive and dynamic maintenance policy
that resembles a CBM policy.
Our idea is described here first in general terms,
and then illustrated along two examplary timelines in Section~\ref{sec:policy-example}.
Afterwards, we give a detailed description of the policy in Section~\ref{sec:operationalprocedure},
and expand on the involved computations in Sections~\ref{sec:costrate} and \ref{sec:optim}.

We denote by $\tnow$ the moment in time that the system is observed,
measured from the moment the present operational cycle has started,
which coincides with starting to use a new system.
Instead of determining a threshold for a degradation signal like in usual CBM approaches,
we determine, for $\tnow$, the adaptive and dynamic time span $\tausnow$
which gives the time from $\tnow$ until the next optimal moment to re-evaluate the system.
%Close to system startup, $\tausnow$ will be large;
%but as components in the system age and some indeed fail over time,
%$\tausnow$ will typically decrease.
A preventive maintenance action is initiated
when $\tnow + \tausnow$ would be earlier than the next planned evaluation time point.
We describe this dynamic and adaptive maintenance policy in more detail in Section~\ref{sec:operationalprocedure}.
%a flow chart for the policy is given in Figure~\ref{fig:procedure}.
%
The optimal time to system re-evaluation $\tausnow$ is determined based on $\Rsysnow(t)$
by minimizing $\gnow(\tau)$, the expected cost rate for the current operational cycle.
%An operational cycle begins with system start-up,
%and ends when either preventive maintenance is carried out,
%or a failure of the system occurs, with subsequent corrective maintenance.
This criterion, as detailed in Section~\ref{sec:costrate},
%only assumes that $\Rsysnow(t)$ is given,
only assumes that a RLD is given,
and thus is not tied to our choice for the component models.
One could also consider alternative optimality criteria to determine $\tausnow$;
when safety is paramount one could, e.g., determine $\tausnow$ such that
the probability of system failure is at most at some very low pre-determined level.
%
Section~\ref{sec:optim} then focuses on our choice for the component models and the resulting form for $\Rsysnow(t)$,
and shows how $\tausnow$ is numerically determined in this situation.
%Numerical examples comparing results for three exemplary failure histories are given in Section~\ref{sec:examples},
%and a small simulation study is presented in Section~\ref{sec:sim}.


\subsection{Two exemplary timelines}
\label{sec:policy-example}

Figure~\ref{fig:timeline} shows two examplary time lines,
starting each with a new system (operational cycle) at $t=0$.
The timeline on the left leads to preventive maintenance,
the timeline on the right to corrective maintenance.

As described in Section~\ref{intro},
the health of the system is re-evaluated at three types of triggers.
Triggers of type (1) are the evenly spaced planned evaluation time points $m\delta$, where $m \in \naturals_0$.
%giving the time granularity at which maintenance actions can be reasonably initiated.***
%Let $\tthresh$ be a threshold time span chosen for practical reasons.
%used to decide when maintenance is initiated, i.e.,
%a time span chosen to reflect operational practicalities.
The times $m\delta$ are indicated by vertical tick lines on the time axis in Figure~\ref{fig:timeline},
where $\delta = 0.1$.
The moments of failure of a component form the other two types of triggers.
A component failure after which the system still functions is a trigger of type (2),
depicted as a (non-filled) circle on the time axis.
A filled circle indicates a trigger of type (3),
which is a component failure that causes the whole system to fail.

For each trigger of type (1) and (2),
the system health is evaluated by calculating $\tausnow$,
the currently optimal time to re-evaluation, %$\tausnow$ is %calculated,
visualized by the thick horizontal bars in Figure~\ref{fig:timeline}.
%The details of calculating $\tausnow$ are described in detail in Sections~\ref{sec:costrate} and \ref{sec:optim}.***
For a type (3) trigger, the system health evaluation just corresponds to acknowledging that the system has failed.
%/ to the diagnosis/insight*** of system failure.
%and calculation of $\tausnow$ is 
For triggers of type (1) and (2),
the calculated $\tausnow$ is compared to the time to the next planned evaluation time point.
When $\tausnow$ is shorter than this time,
preventive maintenance of the system is initiated.

In the timeline on the left of Figure~\ref{fig:timeline},
$\tausnow$ is re-evaluated at times $\tnow = 0.1, 0.2, 0.24, 0.3, 0.36, 0.4$,
where $0.24$ and $0.36$ correspond to component failure times.
Note that the failure at $\tnow = 0.24$ substantially shortens $\tausfun{0.24}$,
such that $0.24 + \tausfun{0.24} < 0.2 + \tausfun{0.2}$,
i.e., the optimal moment of re-evaluation has moved backwards in time.
This illustrates that new insights drawn from a component failure
(regarding the reliability block diagram as well as the component model)
may have considerable influence on the system RLD
and the $\tausnow$ that is derived from it.
The component failure at $\tnow = 0.36$ leads to a further shortened $\tausfun{0.36}$,
such that $\tausfun{0.36} < 0.1 = \delta$.
However, this does not lead to an immediate decision for preventive maintenance,
since $\tausfun{0.36}$ is compared to the time to the next planned evaluation at $t = 0.4$.
The decision for preventive maintenance is instead made at $\tnow = 0.4$,
as only then it is not economical any more to wait for the next planned evaluation at $t = 0.5$.

\begin{figure}
\centering
\begin{tikzpicture}
%[tausnow/.style={|-|, line width=0.5mm}]
[ tausnow/.style={line width=0.8mm},
  pfeil/.style={-latex', line width=0.5mm, shorten >=1mm}] 
\draw[-latex', very thick] (-0.2, -0.5) -- (5.5, -0.5) node[below] {$t$};
\draw (0, -0.4) -- (0, -0.6) node[below] {$0$};
\foreach \x in {1, 2, 3, 4, 5}
 {\draw (\x, -0.4) -- (\x, -0.6) node[below] {$0.\x$};
  \draw[dotted] (\x, -0.4) -- (\x, \x); }
\draw[thick, decorate, decoration=brace] (1, -1) -- (0, -1) node[midway,below] {$\delta$};
\draw[dashed] (0, 0) -- (5,5);
\draw[tausnow] (0  , 0  ) -- (2.3, 0  ) node[at start, left] {$\tau_*^{(0)}$};
\draw[tausnow] (1  , 1  ) -- (2.9, 1  ) node[at start, left] {$\tau_*^{(0.1)}$};
\draw[tausnow] (2  , 2  ) -- (4.2, 2  ) node[at start, left] {$\tau_*^{(0.2)}$};
\draw (2.4, -0.5) circle (0.8mm) edge[dotted] (2.4, 2.4);
%\node[draw,shape=diamond] (f1) at (2.4, -0.4) {}; % edge[dotted] (2.4, 2.4);
\draw[tausnow] (2.4, 2.4) -- (3.9, 2.4) node[at start, left] {$\tau_*^{(0.24)}$};
\draw[tausnow] (3  , 3  ) -- (4.4, 3  ) node[at start, left] {$\tau_*^{(0.3)}$};
\draw (3.6, -0.5) circle (0.8mm) edge[dotted] (3.6, 3.6);
\draw[tausnow] (3.6, 3.6) -- (4.3, 3.6) node[at start, left] {$\tau_*^{(0.36)}$};
\draw[tausnow] (4  , 4  ) -- (4.7, 4  ) node[at start, left] {$\tau_*^{(0.4)}$};
\node (box1) at (1,4) {\parbox{12ex}{\small preventive\\ replacement\\ at $\tnow = 0.4$\\ as $\tau_*^{(0.4)} < \delta$}};
\draw[pfeil] (box1) edge[out=20, in=120] (4,4);
\begin{scope}[xshift=7cm]
\draw[-latex', very thick] (-0.2, -0.5) -- (3.5, -0.5) node[below] {$t$};
\draw (0, -0.4) -- (0, -0.6) node[below] {$0$};
\foreach \x in {1, 2, 3}
 {\draw (\x, -0.4) -- (\x, -0.6) node[below] {$0.\x$};
  \draw[dotted] (\x, -0.4) -- (\x, \x); }
\draw[thick, decorate, decoration=brace] (1, -1) -- (0, -1) node[midway,below] {$\delta$};
\draw[dashed] (0, 0) -- (3,3);
\draw[tausnow] (0  , 0  ) -- (2.3, 0  ) node[at start, left] {$\tau_*^{(0)}$};
\draw[tausnow] (1  , 1  ) -- (2.9, 1  ) node[at start, left] {$\tau_*^{(0.1)}$};
\draw (1.4, -0.5) circle (0.8mm) edge[dotted] (1.4, 1.4);
\draw[tausnow] (1.4, 1.4) -- (2.8, 1.4) node[at start, left] {$\tau_*^{(0.14)}$};
\draw[tausnow] (2  , 2  ) -- (3.3, 2  ) node[at start, left] {$\tau_*^{(0.2)}$};
\draw[fill] (2.7, -0.5) circle (0.8mm) edge[dotted] (2.7, 2.7);
\node (box2) at (1,4) {\parbox{13ex}{\small corrective\\ replacement\\ at $\tnow = 0.27$\\ due to system failure}};
\draw[pfeil] (box2) edge[out=10, in=90] (2.7,2.7);
\end{scope}
\end{tikzpicture}
\caption{Two exemplary timelines for the suggested maintenance policy.
The horizontal bars depict $\tausnow$,
the time to the next economically optimal evaluation of the system RLD,
for subsequent times $\tnow$.
The times $\tnow$ for which $\tausnow$ is evaluated are
at regular intervals $m \delta, m \in \naturals_0$ (trigger 1),
and at the time of component failures (trigger 2),
which are marked by a circle on the time axis.
Failures that lead to the whole system failing (trigger 3) are marked by a filled circle.}
%The calculation of $\tausnow$ is described in detail in Section~***.}
%\begin{tikzpicture}\draw (0,0) circle (0.8mm);\end{tikzpicture}
\label{fig:timeline}
\end{figure}

The timeline on the right of Figure~\ref{fig:timeline} illustrates an operational cycle ending with corrective maintenance.
Here, the component failure at $\tnow = 0.27$ leads to system failure.
The size of $\tausfun{0.2}$ indicates that the component failure that happens at $\tnow = 0.27$
had been deemed very unlikely at time $\tnow = 0.2$,
but then nevertheless occured at $t = 0.27$.


\subsection{Detailed description of the policy}
\label{sec:operationalprocedure}

Each operational cycle starts with a system startup.
At this point in time, the system is in an (as good as) new state,
i.e., all components in the system are (as good as) new and have not aged.
An operational cycle ends when the system is either maintained preventively, %at time $\tnow + \tausnow$ (on an absolute time scale),
or the system fails, which triggers a corrective maintenance action.

We assume that both preventive and corrective maintenance result in an (as good as) new system state by replacing all components,
such that the system begins a new operational cycle always under the same conditions.
At first sight, replacing a complete system when only some parts of it are defective may seem strange.
However, there are many real life examples, including from daily live.
So called complete clustering \citep{2017:oldekeizer}
may be optimal in case of a very strong economic dependence between the different parts of a system,
e.g., due to a very high fixed set-up cost independent of the maintenance activities to be done,
as is the case of installations used in the generation and distribution of oil and gas.
Note that often for safety reasons, the latter installations are monitored continuously.
There may be also technical dependencies between parts of a system that may require combined replacement:
replacing the cassette of a bike often also requires to replace the chain,
and the tires of an airplane are required to have the same thickness.
But also complete systems are replaced because replacing one or more of its parts may take much longer,
which may result in big production losses.
Another reason can be that there may be a high probability that when individual parts are replaced,
other parts will be damaged, requiring their replacement as well.

Replacing not all components during maintenance
would require to account for the different ages of components in the system
at the start of a new cycle.
While this is in principle possible in our approach,
exploring such selective replacement schemes is out of scope for this paper.
Optimizing over replacement schemes is, in our view, a higly interesting topic that merits an investigation of its own.
We will return to this in the conclusion (Section~\ref{sec:outlook}).

Figure~\ref{fig:procedure} depicts the flow chart for the policy,
indicating the decisions that are repeatedly encountered in a single operational cycle.
The flow chart also contains the loop back to the start of a new operational cycle.

%***Consider a set of evaluation time points $\tnow$.
Consider a grid of planned evaluation time points $\{m \delta : m \in \naturals_0\}$,
%where $m = 0, 1, 2, \ldots$,
%and $\delta$ sufficiently small to be negligible on the system lifetime scale.
with $\delta$ a sufficiently small time increment.
During an operational cycle,
$\tausnow$ is usually recalculated at times $\tnow = m \delta$, $m = 0, 1, 2, \ldots$.
%where the time increment $\delta$ should be small enough to be negligible on the scale of system lifetime.
However, if a component $i$ of type $k$ fails before the next $\tnow$ on the grid,
$\tnow$ is set instead to the exact failure time $t_{k,i}$.
Next, we check whether the system has failed due to this component failure.
If the system has failed, a corrective maintenance action is carried out, incurring cost $c_c$.
If the system has not failed, we calculate $\tausnow$ for the current off-grid time $\tnow = t_{k,i}$,
and move back onto the grid afterwards,
unless there is another component failure before the next planned evaluation time.

The calculation of $\tausnow$ for a fixed $\tnow$ resembles the calculation
of the maintenance interval in an age based policy by trading off the costs of preventive and corrective maintenance
(see Section~\ref{sec:costrate} below).
However, because $\tausnow$ is derived from the current $\Rsysnow(t)$,
our policy is adaptive and dynamic like a CBM policy. 
When computation of the RLD is fast,
frequent recalculation of $\tausnow$ is feasible,
%Then, by evaluating the unit cost rate function $\gnow(\tau)$ (see Section~\ref{sec:costrate} below) on a discrete grid,
%$\tausnow$ can be approximated with a simple grid search.
%because we use a fast discrete approximation of the unit cost rate function $\gnow(\tau)$,
%and determine $\tausnow$ with a simple grid search.
%Updating $\tausnow$ frequently is useful because,
and it is useful because,
even in absence of failures of system components,
all components age, and also non-failure adds information to the component models,
reducing uncertainty.
%Both aspects influence $\Rsysnow(t)$ and thus $\tausnow$,
%as is visible in the top panel of Figure~\ref{fig:tauhist1fig2T},
%which gives a numerical example described in more detail in Section~\ref{sec:examples}.
%An alternative would be to update $\tausnow$ only when a component fails,
%but as a (near) real-time update of $\tausnow$ is possible and useful,
%as non-failed components still can be used to update component models ($\nknow$ and $\yknow$) via right-censored observation $\tpnow$.
%
%As can happen for an age-based policy, the derived $\gnow(\tau)$ may also be monotonely decreasing.
%In this case, it is, at the moment $\tnow$, optimal to correctively maintain the system.
%This means that our approach is not only capable of finding the optimal moment for preventive maintenance,
%but also capable of indicating when corrective maintenance is the better option.
%i.e., capable of choosing the maintenance policy.

Once $\tausnow$ is calculated, it is compared to the time span from the current $\tnow$ to the next planned evaluation time.
If the current $\tnow$ is a time on the grid, then $\tnow = m \delta$ for some $m \in \naturals_0$.
The next planned evaluation is then at $(m + 1) \delta$,
and we compare $\tausnow$ to $\delta$.
If $\tausnow < \delta$, a preventive maintenance action is triggered, incurring cost $c_p$,
otherwise we attempt to move to the next planned evaluation point, by setting $m =  (m + 1)$.

If the current $\tnow$ is instead an off-grid time corresponding to failure of component $i$ of type $k$ at time $t_{k,i}$,
after which the system still functions,
then the next planned evaluation time is the next time on the grid.
Let $m$ be such that $m \delta < t_{k,i} < (m + 1) \delta$,
so the next time on the grid is $(m + 1) \delta$,
and we compare $\tausnow$ with $(m + 1) \delta - t_{k,i}$.
If $\tausnow < (m + 1) \delta - t_{k,i}$, a preventive maintenance action is triggered, incurring cost $c_p$,
otherwise we attempt to move to the next planned evaluation point, which is still $(m + 1) \delta$.

In both cases, we decide to maintain the system preventively
when the optimal time to re-evaluation is closer to now than the next time we intend to re-evaluate the system RLD, %$\tausnow$,
since otherwise, in waiting to the next planned evaluation time,
one would exceed the optimal system failure risk 
as determined by the trade-off made via $\gnow(\cdot)$. %***drop $\gnow(\cdot)$ here, like: ...trade-off that results in $\tausnow$.
When $\tausnow$ is instead further in the future than the next planned evaluation time,
we can safely continue to the next planned evaluation time,
without risking to miss the cost-optimal moment to preventively maintain the system.

\begin{figure}
\centering
\begin{tikzpicture}
[node distance=6mm,
 font=\footnotesize,
 box/.style={draw, rectangle, thick, minimum height=6mm, minimum width=15mm},
 box2/.style={draw, rectangle, thick, minimum height=6mm, minimum width=28mm, rounded corners=2mm},
 dec/.style={draw, shape=diamond, thick, minimum height=4mm, minimum width=4mm},
 pfeil/.style={-latex', very thick}]
\node (startup) [box] {\parbox{18ex}{\centering system start-up\\ $\tnow = 0, m= 0$\\ initial params.\\ $\nkz, \ykz$}};
\node (mincrem) [box, below=57mm of startup] {$m = m + 1$};
\node (cptfail) [dec, below right=9mm and 17mm] {\parbox{10ex}{\centering $\exists\, t_{k,i}$ s.t.\\ \rule{0ex}{2ex} \\ ?}};
%$(m-1)\delta$\\ $ < t_{k,i} < $\\ $m\delta$}};
\node at (cptfail) {$\scriptstyle \tnow < t_{k,i} < (m+1)\delta$};
%\coordinate (mcoord) at ($(startup) * (0, -2)$) {};
\node (nfltnow) [box, below=of cptfail] {\parbox{21ex}{\centering $\tnow = (m + 1)\delta$,\\ update parameters\\ to $\nknow, \yknow$}};
\node (yfltnow) [box, right=of cptfail] {\parbox{21ex}{\centering $\tnow = t_{k,i}$,\\ update parameters\\ to $\nknow, \yknow$}};
\node (sysfail) [dec, below=5.6mm of yfltnow] {\parbox{8ex}{\centering system\\ failed?}};
\node (tausmde) [box, below=of nfltnow] {\parbox{21ex}{\centering determine $\tausnow$}};
\node (taustki) [box, below=of sysfail] {\parbox{21ex}{\centering determine $\tausnow$}};
\node (tausdel) [dec, below=of tausmde] {\parbox{10ex}{\centering $\tausnow$\\ \rule{2ex}{0ex}$< \delta$?}};
\node (taudtki) [dec, below=of taustki] {\parbox{10ex}{$\tausnow <$\\ \phantom{?}}};
\node at ($(taudtki) + (0,-2mm)$) {$\scriptstyle (m + 1) \delta - t_{k,i}$?};
\node (corrmnt) [box, right=8.5mm of sysfail] {\parbox{21ex}{\centering corrective repair\\ at $\tnow$ (cost $c_c$)}};
\node (prevmnt) [box, right=of taudtki] {\parbox{21ex}{\centering preventive repair\\ at $\tnow$ (cost $c_p$)}};
\node (updatep) [box, below=7.5mm of corrmnt] {\parbox{21ex}{\centering set initial params.\\ for next cycle as\\ $\nkz = \nknow$\\ $\ykz = \yknow$}};
%\coordinate (nocoord) at ($(taudtki) + (-1.8,-1.3)$) {};
\coordinate (nocoord) at ($(taudtki) + (-1.8, 1.6)$) {};
\coordinate (yscoord) at ($(tausdel) + ( 0  ,-1.5)$) {};
\coordinate (nxtcoor) at ($(updatep) + ( 1.8, 0  )$) {};
\draw[pfeil] (startup) |- ($(cptfail.west) + (0.1, 0.1)$);
\draw[pfeil] (mincrem) |- ($(cptfail.west) + (0.1,-0.1)$);
\draw[pfeil] (cptfail) -- (nfltnow) node[midway, left] {no};
\draw[pfeil] (cptfail) -- (yfltnow) node[midway, above] {yes};
\draw[pfeil] (yfltnow) -- (sysfail);
\draw[pfeil] (nfltnow) -- (tausmde);
\draw[pfeil] (sysfail) -- (taustki) node[midway, left] {no};
\draw[pfeil] (sysfail) -- (corrmnt) node[midway, above] {yes};
\draw[pfeil] (tausmde) -- (tausdel);
\draw[pfeil] (taustki) -- (taudtki);
\draw[pfeil] (tausdel) -| (mincrem) node[at start, above left] {no};
\draw[pfeil] (taudtki) -| (nocoord) node[near start, above] {no} -- ($(mincrem) + (0,7mm)$); %($(0,1cm)*(nocoord) + (1cm,0)*(mincrem)$); %
\draw[pfeil] (taudtki) -- (prevmnt) node[midway, above] {yes};
\draw[pfeil] (tausdel) -- (yscoord) -| (prevmnt) node[at start, above right] {yes};
\draw[pfeil] (prevmnt) -- (updatep);
\draw[pfeil] (corrmnt) -- (updatep);
\draw[pfeil, dotted] (updatep) -- (nxtcoor) |- (startup) node[midway, above left] {start new cycle};
%
%\node (startup) [box]                       {\parbox{20ex}{\centering system start-up\\ (all comp.\ new)}};
%\node (sysfail) [box, right=of startup]     {system failed?};
%\node (tauchck) [box, below=6mm of sysfail] {$\tausnow < \delta$?};
%\node (gotonxt) [box, below=6mm of tauchck] {\parbox{17ex}{\centering continue to\\ $\tnow = (m+1)\delta$}};
%\draw [thick, dashed, rounded corners=2mm] ($(gotonxt.south east) + (2mm,-1mm)$) rectangle ($(sysfail.north west) + (-5mm,1mm)$) node[above left=0mm and -33mm] {each $\tnow = m \delta, m \in \naturals_0$};
%\node [above left=0mm and -8mm of startup] {$t=0$};
%\node (correct) [box, right=of sysfail, minimum width=18mm] {\parbox{14ex}{\centering corrective\\ maintenance\\ (cost $c_c$)}};
%\node (prevent) [box, right=of tauchck, minimum width=18mm] {\parbox{14ex}{\centering preventive\\ maintenance\\ (cost $c_p$)}};
%\node (update) [box, below right=10mm and -4mm of prevent, minimum width=10mm] {\parbox{12ex}{\centering update parameters}};
%\draw[pfeil] (startup) -- (sysfail);
%\draw[pfeil] (sysfail) -- (tauchck) node[midway, right] {no};
%\draw[pfeil] (tauchck) -- (gotonxt) node[midway, right] {no};
%\draw[very thick] (gotonxt.west) -| ($(sysfail.west) + (-3mm,0mm)$);
%\draw[pfeil] (sysfail) -- (correct) node[midway, above] {yes};
%\draw[pfeil] (tauchck) -- (prevent) node[midway, above] {yes};
%\draw[pfeil] (correct) -| (update);
%\draw[pfeil] (prevent) -| (update);
%\draw[pfeil, dotted] (update) -| (startup) node[midway, below right] {start new cycle};
\end{tikzpicture}
\caption{Flowchart for the proposed maintenance policy,
showing a single operational cycle,
including the loop back to system startup that starts a new cycle as a dashed arrow.}
%For simplicity of presentation, we exclude the possibility that several component failure times
%happen between the evaluation times $(m-1) \delta$ and $m \delta$.}
\label{fig:procedure}
\end{figure}

At the end of each operational cycle, i.e.,
after maintaining the system either preventively or correctively,
the component information gained during the completed operational cycle is added to the component knowledge base
by updating the parameters $\nk$ and $\yk$ accordingly,
%using failure and censoring times as of the end of the cycle.
such that $\nkz$ and $\ykz$ for the next operational cycle
account for the (censored) lifetimes observed during the previous operational cycles.
%We assume that both planned and unplanned maintenance action takes negligible time,

%
%As we use a one-cycle criterion, the time both preventive and corrective maintenance actions require
%does not matter directly for the determination of $\tausnow$.
%Longer downtime and other negative effects of corrective maintenance
%are instead accounted for through the costs parameters $c_p$ and $c_c$ used in the trade-off leading to $\tausnow$,
%as described next.


\subsection{Expected operational cycle cost rate}
\label{sec:costrate}

To determine the optimal time to maintenance $\tausnow$,
we use the expected one-cycle cost rate 
\citep{1984:ansell-bendell-humble,1996:mazzuchi-soyer,2006:coolen-schrijner-coolen}
as the unit cost rate.
As discussed in the literature review (Section~\ref{sec:literature}),
we cannot use the renewal theory based way to calculate the unit cost rate as often encountered
in CBM literature \citep[e.g.,][]{2013:si-et-al,2011:kim-et-al}.
%we use the one-cycle cost rate as it is the most appropriate for our situation,
%where $\tausnow$ changes frequently.

Let $c_p$ be the cost of preventive maintenance, and $c_c$ the cost of corrective (unplanned, breakdown) maintenance, where $c_p < c_c$.
Usually, set-up costs, the actual maintenance costs (labor, new components) and downtime costs are much higher for corrective maintenance,
such that often, $c_p \ll c_c$.

Let $\Tsysnow$ be the random variable corresponding to $\Rsysnow(t)$,
giving the (random) system failure time for times $t > \tnow$,
conditioned on all information available at $\tnow$,
and $\fsysnow(t)$ be the corresponding probability density.
Note that due to the conditional nature of $\Rsysnow(t)$,
which gives the probability that the system functions at time $t$ given that it functions at $\tnow$,
also $\fsysnow(t)$ is a conditional probability density,
with support $\{t \in \reals : t > \tnow \}$.
%In particular, we have $\fsysnow(t) = 0$ for all $t \le \tnow$.

Conditional on the realization $\tsysnow$ of $\Tsysnow$, the unit cost rate 
for planning maintenance at $\tau$ time units after $\tnow$ is 
\begin{align}
g(\tau \mid \Tsysnow = \tsysnow) &=
\begin{cases}
%c_p / \tau               & \text{if } \tsysnow \ge \tnow + \tau \\
%c_u / (\tsysnow - \tnow) & \text{if } \tnow < \tsysnow  <  \tnow + \tau \,.
c_p / (\tnow + \tau) & \text{if } \hspace{7.4ex} \tsysnow \ge \tnow + \tau \\
c_c /  \tsysnow      & \text{if }        \tnow < \tsysnow  <  \tnow + \tau \,.
\end{cases}
\label{eq:gtau}
\end{align}
Note that here $\tau$ is a time span starting from $\tnow$,
and so is a time on a prospective time scale, where $0$ corresponds to $\tnow$.
In contrast, $\tnow$ and $\tsysnow$ are times on an absolute time scale,
where $0$ corresponds to the time of a new system start-up.
Taking the expectation over $\Tsysnow$ in \eqref{eq:gtau} leads to the expected operational cycle cost rate
\begin{align}
\gnow(\tau) &= \E[g(\tau \mid \Tsysnow)] \nonumber\\ 
            &= \frac{c_p}{\tnow + \tau} \Rsysnow(\tnow + \tau)
              + c_c \int_0^\tau \frac{1}{\tnow + \theta} \fsysnow(\tnow + \theta) \dd \theta\,.
\label{eq:gtnowtau}
\end{align}
Here, we assume that the time required to execute preventive and corrective maintenance
can be neglected when compared with the duration of one operational cycle.
%
$\gnow(\tau)$ gives the cost rate that we expect to incur over the complete operational cycle
when, at time $\tnow$, we schedule preventive maintenance for time $\tnow + \tau$,
so $\gnow(\tau)$ trades off the risk of system failure before $\tnow + \tau$
with the gains we could reap when the system survives until $\tnow + \tau$.
Figure~\ref{fig:ghist1fig2T} shows $\Rsysnow(t)$ and the corresponding $\gnow(\tau)$ for a number of current times $\tnow$
in a numeric example which is described in more detail in Section~\ref{sec:examples}.

\begin{figure}
\includegraphics[width=\textwidth]{ghist1fig2T}
\caption{System reliability and unit cost rate functions for $\tnow = 0,2,4,6,8$
for the system depicted in Figure~\ref{fig:brakesys-layout}
and prior component models as given in Table~\ref{tab:priorparams},
assuming the failure history $P_2 = 3$, $P_3 = 4$, $C_2 = 6$, $C_3 = 7$, $H = 8$,
and costs $c_u = 1$, $c_p = 0.2$.}
\label{fig:ghist1fig2T}
\end{figure}

Note that for the integral in \eqref{eq:gtnowtau} to exist,
%$\E[1/\Tsysnow \mid \Tsysnow > \tnow]$ must exist,
$\E[1/(\Tsysnow - \tnow)]$ must exist,
which is a condition for $\fsysnow(\tnow + \theta)$ for values of $\theta$ close to $0$ \citep{2006:coolen-schrijner-coolen}.
%This restriction probably plays a role in the CBM literature's reluctance to adopt the one-cycle criterion.
For our approach, however, this condition is not relevant,
as we use a discrete approximation of $\Rsysnow(t)$ and $\fsysnow(t)$, respectively
(see Section~\ref{sec:optim} below).

$\tausnow := \arg\min \gnow(\tau)$
is the cost-optimal time to do maintenance as of $\tnow$,
incorporating all that is known at current time $\tnow$.
%The lower index $\tnow$ of $\tausnow$ emphasizes that
%depends on $\tnow$, and also on the history of the system.
%illustrates the derivation of $\tausnow$ .
%Minimizing $\gnow(\tau)$ leads to $\tausnow$,
%the optimal moment of maintenance on a prospective time scale.
The minimal expected operational cycle cost rate $\gstarnow := \gnow(\tausnow)$ corresponding to $\tausnow$
thus is, as of $\tnow$, the lowest cost rate for the current operational cycle
that we can expect to attain.
%attainable for the span of the remaining lifetime of the system,
%i.e., for what happens beyond current time $\tnow$.
%$\gstarnow$ increases over time since 
%***Also of interest is the total expected cost rate as of $\tnow$, denoted by $\gtotalnow$.
%This gives the expected unit cost rate calculated over the full time since system start-up,
%under the assumption that maintenance is scheduled for $\tstarnow := \tnow + \tausnow$.
$\gstarnow$ generally decreases with $\tnow$ since for larger $\tnow$, maintenance costs are spread over a longer time period.
However, when many (or few but critical) components fail, $\gstarnow$ may increase,
since then, system failure becomes much more likely, leading to a higher risk of incurring the corrective maintenance cost $c_c$.
Figure~\ref{fig:tauhist1fig2T} illustrates these relations for the scenario depicted in Figure~\ref{fig:ghist1fig2T}.
%
As can be seen there and in the numerical examples in Section~\ref{sec:examples},
for continuous component lifetime distributions,
$\tausnow$ is a smooth function of $\tnow$ during periods where no component fails,
with discrete jumps at component failure times.
The size of the jump depends on the relevance of the failed component
for the reliability of the system at the moment that the component failure occurs.

Note that common cause failures, i.e., failures of several components at the same time due to a shared root cause,
are currently not accounted for in our model. %, as we assume component failures to be mutually independent
%given the component model parameters.
Such simultaneous failures will seriously influence the RLD and thus $\tausnow$.
We will return to this in the conclusion (Section~\ref{sec:outlook}).
%
%Next, we describe how to numerically evaluate \eqref{eq:gtnowtau} and determine $\tausnow$
%for our choice of component model.

\begin{figure}
%\centering
\includegraphics[width=\textwidth]{tauhist1fig2T}
\caption[t]{$\tausnow$, $\tstarnow = \tnow + \tausnow$ and $\gstarnow = \gnow(\tausnow)$
%and the total unit cost rate $\gtotalnow$ as functions of $\tnow$
for the system depicted in Figure~\ref{fig:brakesys-layout},
assuming the failure history $P_2 = 3$, $P_3 = 4$, $C_2 = 6$, $C_3 = 7$, $H = 8$,
and costs $c_u = 1$, $c_p = 0.2$.}
\label{fig:tauhist1fig2T}
\end{figure}


\subsection{Numerical optimization of the expected operational cycle cost rate}
\label{sec:optim}

For our choice of the component models,
the resulting $\Rsysnow(t)$ as of \eqref{eq:sysrel-tnow}
is closed-form, but not tractable analytically.
Therefore,
we evaluate $\Rsysnow(t)$ on a dense grid of time points $t$
covering the time span from $\tnow$ to some horizon time $\tnow + t_h$.
%(The choice of $t_h$ is discussed in more detail below.)
We calculate $\gnow(\tau)$ on the same grid,
and determine $\tausnow$ by minimizing over this grid.
Here, we approximate the integral in \eqref{eq:gtnowtau} numerically by a sum,
where $\fsysnow(t)$ is calculated via differences of $\Rsysnow(t)$.
%$\gnow(\tau)$ is likewise evaluated on a dense grid,
%and $\tausnow = \arg\min \gnow(\tau)$ is determined as the minimum over the grid.

The horizon time span $t_h$ can be determined by practical considerations,
but generally should reach far enough into the tail of $\Rsysnow(t)$
such that $\gnow(\tau)$ is not cut off before it reaches its minimum. %one can be sure to find $\tausnow$.
Note that according to \eqref{eq:gtnowtau},
the calculation of $\gnow(\tau)$ requires
evaluations of $\Rsysnow(t)$ and $\fsysnow(t)$ only for $\tnow < t \le \tnow + \tau$.
The choice of $t_h$ is thus not critical,
in the sense that $\gnow(\tau)$ evaluations will be correct even if $t_h < \tausnow$.

We have implemented our method in \textsf{R} \citep{R},
using the package \texttt{ReliabilityTheory} \citep{2016:aslett-RT}
to calculate survival signatures
(code available on request).
As the parameter update step is closed-form,
recalculation of $\tausnow$ is a matter of seconds,
allowing for a near real-time update.
%
%***calculations in \textsf{R} can be vectorized and are thus fast.
%***if min $\gnow(\tau)$ at $\tnow + t_h$, corrective maintenance is optimal: policy switch!
%***if $\tau_*^{(\tnow)}$ larger than certain horizon, than conclude corrective policy


\section{Further elaboration of the case study}
\label{sec:examples}

We now explain in more detail the case study used for the previous illustrative figures
(Figures~\ref{fig:ghist1fig2T} and \ref{fig:tauhist1fig2T}).
We use the simplified automotive brake system as introduced in Section~\ref{sec:sysrel},
see Figure~\ref{fig:brakesys-layout} for the reliability block diagram.

The cost parameters are chosen as $c_c = 1$, $c_p = 0.2$, % in all examples,
i.e., corrective maintenance is five times more costly than preventive maintenance.
This makes it worthwile to aim for a preventive repair of the system,
but corrective maintenance will not be avoided at all cost.
%Indeed, $c_u / c_p = 5$ means that risking corrective maintenance is worthwile
%if one can expect to run the system five times longer
%than the preventive maintenance interval.

The choice of prior parameters for the four component models are discussed in Section~\ref{sec:ex-prior}.
Then, our maintenance policy is illustrated under three different exemplary failure time histories,
described in Sections~\ref{sec:ex-1}, \ref{sec:ex-2} and \ref{sec:ex-3}.

\begin{table}
\centering
\begin{tabular}{crrrrrrr}
  \toprule
$k$ & $\beta_k$ & $\E[T_i^k]$ & $\ykz$ & $\nkz$ \\
  \midrule
% test prior
%M & $2.5$ & $5$\rule{1.5ex}{0ex} & $75.4$ & $2$\rule{1ex}{0ex} \\
%H & $1.2$ & $2$\rule{1.5ex}{0ex} & $ 2.5$ & $1$\rule{1ex}{0ex} \\
%C & $2  $ & $8$\rule{1.5ex}{0ex} & $81.5$ & $1$\rule{1ex}{0ex} \\
%P & $1.5$ & $3$\rule{1.5ex}{0ex} & $ 6.1$ & $1$\rule{1ex}{0ex} \\
% prior 1
C & $2  $ & $ 3$\rule{1.5ex}{0ex} & $ 11.5$ & $1$\rule{1ex}{0ex} \\
H & $1  $ & $10$\rule{1.5ex}{0ex} & $ 10.0$ & $1$\rule{1ex}{0ex} \\
M & $2.5$ & $ 8$\rule{1.5ex}{0ex} & $244.1$ & $1$\rule{1ex}{0ex} \\
P & $1.5$ & $ 5$\rule{1.5ex}{0ex} & $ 13.0$ & $1$\rule{1ex}{0ex} \\
  \bottomrule
\end{tabular}
\caption{Prior parameters for the four component types as used in the numerical examples.}
%The corresponding expected failure behaviour is visualized in Figure~\ref{fig:compprior1fig1}.}
\label{tab:priorparams}
\end{table}


\subsection{Prior component models}
\label{sec:ex-prior}

The prior component model parameters chosen for this case study
are given in Table~\ref{tab:priorparams}.
Note that the scale parameter for component type H is $1$,
such that the failure times for this component type are assumed to follow an Exponential distribution.
In cases where only $H$ and one or both of $P_3$ and $P_4$ are functioning,
the system reliability function will thus be close to exponential,
making preventive maintenance unattractive.

Note also that $\nkz = 1$ for all four component types,
indicating very weak prior knowledge.
Observed failure times thus have a strong influence on the updating of $\ykz$,
such that component models are primarily data-driven.
This offers the advantage of quick adaption
for the case that an expert's assessment of the component MTTF turns out to be inadequate.
%things could look different for stronger expert information!
When $\nkz$ and $\ykz$ also include test data,
then $\nkz$ will typically be larger,
since the update of $n_k$ involves the number of observed failures
(see Section~\ref{sec:weibull}).

The expected component failure behaviour obtained by the parameter choices in Table~\ref{tab:priorparams}
is visualized in Figure~\ref{fig:compprior1fig1},
showing the prior predictive reliability function for a single component of each type.
Note that these functions are not plain Weibull reliability functions,
but weighted averages over Weibull reliability functions according to the prior distribution over the scale parameter $\lambda_k$.
Judging the adequacy of the prior parameter choices thus requires the sampling model to be correct.
Methods for checking the sampling model are discussed, e.g., in \citet{2000:bayarri-berger}.
Figure~\ref{fig:compprior1fig1} indicates that, e.g., the median prior expected failure times are
$2.2$, $4.1$, $6.3$, and $3.1$ for type $C$, $H$, $M$, and $P$ components, respectively.
While it is a priori quite unlikely that any type $C$ component will survive past time $6$,
we expect 25\% of type $H$ components to survive past time $10$.

\begin{figure}
\centering
\includegraphics[width=\textwidth]{compprior1fig1}
\caption{Prior predictive reliability functions corresponding to the values from Table~\ref{tab:priorparams},
%for the four component types of the system depicted in Figure~\ref{fig:brakesys-layout},
showing the expected failure behaviour.}
\label{fig:compprior1fig1}
\end{figure}


\subsection{Failure history 1: failures later than expected}
\label{sec:ex-1}

In the first failure history, we observe
$P_2 = 3$, $P_3 = 4$, $C_2 = 6$, $C_3 = 7$ and $H = 8$,
and in following the system until $\tnow = 8$,
all other failure times are right-censored.
This is the failure history already used for Figures~\ref{fig:ghist1fig2T} and \ref{fig:tauhist1fig2T}.
Here, failures observed for type $C$ are thus considerably later than expected,
while for type $P$, the two observed failures are close to the assumed MTTF. 
Both failure at time $8$ for $H$ and non-failure by time $8$ for $M$ are nothing unusual.
In total, due to the late type $C$ failures, prior assumptions turn out to be somewhat pessimistic.

This is reflected in Figure~\ref{fig:tauhist1fig2T},
where $\tausnow$ (left panel) first decreases quickly but then tends to decrease more slowly during failure-free periods,
which is a consequence of adjusting prior expectations for type $C$ components
due to their unexpected reliability.
Large drops of $\tausnow$ happen only at the failures of components $P_2$, $P_3$ and $H$
at times 3, 4 and 8, respectively,
as $\Rsysnow(t)$ changes considerably at these moments.
$\tausnow$ falls below $\delta = 0.1$ at time $\tnow = 6$,
triggering preventive maintenance of the system.
$\tausnow = 1.60$ at $\tnow = 0$, this corresponds to the optimal maintenance time
in an age-based policy that takes the system as monolithic. 
Since $\tausnow$ does not change dramatically most of the time, $\tstarnow$,
the cost-optimal moment of maintenance on the absolute timescale (center panel),
steadily increases, with small drops corresponding to drops in $\tausnow$,
showing that as time progresses, accumulated information indicates that maintenance can be postponed.
Consequently, $\gstarnow$ (right panel) is steadily decreasing,
with small upward jumps at the failures of $P_2$ and $P_3$.
%so the drops in $\tausnow$ do not entirely outweigh the higher likelihood of system failure
%due to the failure of these components.
%\textbf{***show $\Rsysnow(\tausnow)$???}
At $\tnow = 6$, the expected operational cycle cost rate is $0.033$,
and so considerably lower than the operational cycle cost rate of $0.2 / 1.6 = 0.125$
that one obtains by maintaining the system preventively at time 1.6,
as suggested by the age-based policy that disregards the information from component monitoring.
In this example, using the component failure information allows to lower the cost rate by 73\%.

\begin{figure}
\includegraphics[width=\textwidth]{tauhist1fig3T}
\caption{Illustration of the gains due to the parameter update at each $\tnow$,
assuming the failure history $P_2 = 3$, $P_3 = 4$, $C_2 = 6$, $C_3 = 7$, $H = 8$.}
\label{fig:tauhist1fig3T}
\end{figure}

The effect of the continuous parameter update in our policy is illustrated in Figure~\ref{fig:tauhist1fig3T}.
There, the results as seen already in Figure~\ref{fig:tauhist1fig2T} are compared
to the case when component parameters are not updated during the operation of the system,
while still accounting for the fact that components age,
and that a failed component changes the system layout.
%***so dynamic but not adaptive?***
We see that $\tausnow$ (left panel) is decreasing more quickly in the no update case,
with the jumps at component failure times being less pronounced than in the update case.
Without the continuous update, residual system lifetime continues to be underestimated,
and maintenance is scheduled too early,
%$\tausnow$ (left panel) is consistently lower for the non-adaptive algorithm,
%same with $\tstarnow$ (center panel).
%This is because the adaptive algorithm recognizes that system components are more reliable than expected,
leading to higher expected operational cycle cost rates $\gstarnow$ (right panel).
$\tausnow$ for the policy without parameter update drops below $\delta = 0.1$ at $\tnow = 4$,
while maintenance is triggered for the algorithm with parameter update at $\tnow = 6$.
The parameter update thus allows to run the system for two more time units,
corresponding to a unit cost rate of $0.2 / 6 = 0.033$ instead of $0.2 / 4 = 0.05$,
a reduction of 33\%.


\subsection{Failure history 2: failure times as expected}
\label{sec:ex-2}

In a second failure history, we observe the failure times $C_2 = 1$, $C_2 = 3$, $P_2 = 0.5$, $P_3 = 1.5$,
and follow the system until $\tnow = 2.5$.
Here, components of type $C$ behave more or less as expected,
while type $P$ components fail slightly earlier than expected.
Overall, prior assumptions are reasonably in line with observed failure behaviour.
%with a slight tendency to earlier failures.
Consequently, the difference between the continuous update case and the no update case is small,
as visualized in Figure~\ref{fig:tauhist2fig3T}.
There, $\tausnow$ values for the policy with parameter updates vary around those obtained for the policy without parameter updates,
and a similar behavior for $\tstarnow$ and $\gstarnow$ is observed.
With $\delta = 0.1$, preventive maintenance is carried out at times $2.6$ and $2.5$ for the two policies, respectively.
This shows that for this failure history, in which there is not much to be gained by the parameter update,
$\tausnow$ behaviour is dominated by the update of the system reliability block diagram.
Here, any differences due to the update are very small in comparison to the gains offered by the continuous update
in case of the first failure history
(note the differing vertical and horizontal scales in Figures~\ref{fig:tauhist1fig3T} and \ref{fig:tauhist2fig3T}).

\begin{figure}
\includegraphics[width=\textwidth]{tauhist2fig3T}
\caption{Effect of the continuous parameter update for the second failure history ($C_2 = 1$, $C_2 = 3$, $P_2 = 0.5$, $P_3 = 1.5$)
which is more or less in line with expectations.
%containing slightly early failures of type $P$ components.
With observed failure behaviour more or less in line with expectations,
both the policies with and without continuous update perform similarly,
unlike in case of the first failure history as depicted in Figure~\ref{fig:tauhist1fig3T},
where updating gave considerable gains.}
\label{fig:tauhist2fig3T}
\end{figure}


\subsection{Failure history 3: failures earlier than expected}
\label{sec:ex-3}

%The third failure history contains very late failures for all component types
%by assuming $C_2 = 10$, $C_3 = 12$, $H = 14$, $P_2 = 8$, $P_3 = 9$,
%and following the system until $\tnow = 14$.
%The late failure times for type $C$ components are especially surprising,
%with the prior predictive reliability being $0.01$ at time $10$
%(see top left panel of Figure~\ref{fig:compprior1fig1}).
%Here, the situation is more extreme than in the first failure history.
%Figure~\ref{fig:tauhist3fig3} shows that, after a first short dip,
%$\tausnow$ for the continuous update model steadily increases until the first observed failure at $\tnow = 8$.
%For the model without update, $\tausnow$ steadily decreases instead,
%ignoring any information on component behaviour to be gained during the current system run.
%Assuming $\tthresh = 0.5$, the method without parameter update would recommend to maintain the system preventively at $\tnow = 8.4$,
%with only $P_2$ having failed by then, and all other components up and running,
%corresponding to a considerable waste of component lifetime.
%
%\begin{figure}
%\includegraphics[width=\textwidth]{tauhist3fig3}
%\caption{Effect of the continuous parameter update for the third failure history,
%containing very late failures for all component types ($C_2 = 10$, $C_3 = 12$, $H = 14$, $P_2 = 8$, $P_3 = 9$).}
%\label{fig:tauhist3fig3}
%\end{figure}
%
The third failure history contains very early failures
with $C_2 = 0.1$, $C_3 = 0.2$, $C_4 = 0.8$ $H = 1$, $P_2 = 0.3$, $P_3 = 0.4$,
and we follow the system until $\tnow = 1$,
see Figure~\ref{fig:tauhist4fig3T}.

The policy without parameter update leads to consistently higher $\tausnow$ values,
while the policy with parameter update takes the surprisingly early failures into account.
Unlike in the other two failure histories, $\tstarnow$ is generally decreasing,
indicating that maintenance should happen much earlier than expected initially.
The policy with parameter update consistently suggests to maintain earlier than the policy without update,
but for both models, $\tausnow$ drops below $\delta = 0.1$ upon the failure of H at $\tnow = 1$,
obtaining an operational cycle cost rate of $0.2$.
In contrast, 
the age-based policy that disregards the information from component monitoring
would have scheduled maintenance for time $1.6$,
but it is very unlikely that the system will still function by that time.
If the system fails between time 1 and time 1.6,
then the age-based policy would lead to an operational cycle cost of $0.625$ to $1$,
much higher than for our policy.
%expected cycle costs are underestimated without update***

Like in failure history 2, the structural information (which components still function?)
seems to deliver a major part of the gains for our CBM policy,
whereas the parameter update seems most effective in situations where
component failures happen later than expected.
These findings are confirmed by results from a small simulation study described next.

\begin{figure}
\includegraphics[width=\textwidth]{tauhist4fig3T}
\caption{Effect of the continuous parameter update for the third failure history,
containing very early failures ($C_2 = 0.1$, $C_3 = 0.2$, $C_4 = 0.8$ $H = 1$, $P_2 = 0.3$, $P_3 = 0.4$).}
\label{fig:tauhist4fig3T}
\end{figure}


\section{A small simulation study}
\label{sec:sim}

%A systematic study based on simulation failure times is subsequently presented in Section~\ref{sec:sim}.
We conducted a small simulation study to illustrate the performance of our policy
as compared to corrective and age-based policies.
To be able to discern the effect of parameter learning over several operational cycles,
we simulated component failure times for five subsequent operational cycles for each simulation repetition.
We thus can compare policies that use information gained in previous operational cycles
%We thus can compare policies that carry information gained in one operational cycle to the next operational cycle
with policies that do not update component parameters at the end of an operational cycle.
We also distinguish condition-based maintenance policies with and without continuous parameter update
as in Section~\ref{sec:examples}.
In total, we compare the six maintenance policies listed in Table~\ref{tab:sim-abbrev},
which also gives the abbreviations used in Figures~\ref{fig:br1sim1fig5Tt01}, \ref{fig:br1sim2fig5Tt01} and \ref{fig:br1sim3fig5Tt01}.

\begin{table}
\centering
\begin{tabular}{ll}
  \toprule
Policy & Description \\
  \midrule
CBM-cpu & Condition-based maintenance (continuous parameter update)\\
CBM-epu & Condition-based maintenance (end of cycle parameter update)\\
CBM-npu & Condition-based maintenance (no parameter update)\\
ABM-epu & Age-based maintenance (end of cycle parameter update)\\
ABM-npu & Age-based maintenance (no parameter update)\\
CM      & Corrective maintenance\\
  \bottomrule
\end{tabular}
\caption{The six maintenance policies compared in the simulation study.}
\label{tab:sim-abbrev}
\end{table}

We compare these policies by three simple performance measures,
calculated for each repetition. %, which contains five operational cycles.
$\esys$ is the number of system failures that occured when applying the policy,
This is an integer between zero and five, as there are five operational cycles per repetition.
$\mrsys$ is the mean realized system runtime,
i.e., the average of the realized operational cycle length,
taken over the five operational cycles per repetition.
%(I called this 'cycle length' before),
$\bar{g}$ is the mean realized operational cycle cost rate obtained over the five operational cycles.
This is calculated as the sum of preventive and corrective maintenance costs incurred in the five operational cycles,
divided by the sum of the five realized operational cycle lengths.
%This is a weighted average of the realized operational cycle cost rate for each of the five operational cycles,
%with weights proportional to the corresponding realized operational cycle length.
%This equals 

We use again the simplified automotive brake system depicted in Figure~\ref{fig:brakesys-layout},
together with the prior parameters for the component models as given in Table~\ref{tab:priorparams},
and assume that $c_c = 1$ and $c_p = 0.2$ like before.
We study three data scenarios A, B and C, corresponding to failure times as expected, as well as
failure times earlier and later than expected.

\begin{figure}
\includegraphics[width=\textwidth]{br1sim1fig5Tt01ac}
\caption{Simulation study with 20 repetitions, where each repetition consists of 5 subsequent operational cycles.
Case A: $\lambda_k = \ykz$.}
\label{fig:br1sim1fig5Tt01}
\end{figure}

\subsection{Case A: failure times as expected, $\lambda_k = \ykz$}
\label{sec:case1}

For this scenario, we simulated component failure times according to prior expectations,
i.e., component failure times were simulated as Weibull failure times
with the scale parameter $\lambda_k$ chosen as equal to the expert's prior guess $\ykz$
for all four component types.
The results are depicted in Figure~\ref{fig:br1sim1fig5Tt01}.

Obviously, the corrective maintenance policy (CM) results in $\esys = 5$ for all 20 repetitions,
as it cannot lead to preventive maintenance.
In contrast, all condition-based (CBM) and age-based policies can avoid system failures entirely (with a single exception for ABM-epu).
Mean system runtimes are longest for CM because it runs the system until failure,
but this leads to the highest mean cost rates.
The condition-based maintenance policy with continuous parameter update (CBM-cpu) obtains the lowest mean cost rates,
since it leads to the longest system runtimes without system failures.
CBM-epu and CBM-npu have slightly higher cost rates,
as they lead to shorter system runtimes.
The age-based policy without parameter update (ABM-npu)
always schedules maintenance for $\tau_*^{(0)} = 1.6$.
As the system does not fail before this time in any of the five operational cycles within the 20 repetitions,
the mean realized system runtime is constantly $1.6$, with mean cost rate $0.2 / 1.6 = 0.125$.

\begin{figure}
\includegraphics[width=\textwidth]{br1sim2fig5Tt01ac}
\caption{Simulation study with 20 repetitions, where each repetition consists of 5 subsequent operational cycles.
Case B: $\lambda_k = 0.5\ykz$.}
\label{fig:br1sim2fig5Tt01}
\end{figure}

\subsection{Case B: failure times earlier than expected, $\lambda_k = 0.5 \ykz$}
\label{sec:case2}

For the second scenario,
component failure times were simulated as being earlier than expected,
by choosing $\lambda_k$ for all component types as half the size of the expert's prior guess,
i.e., $\lambda_k = 0.5 \ykz$ for all $k =$ C, H, M, P.
%($0.5 \times$ MTTF) than expected.
The results are depicted in Figure~\ref{fig:br1sim2fig5Tt01}.

%Here, most condition-based and all age-based policies cannot avoid system failures entirely.
Here, only CBM-cpu can avoid system failures entirely.
Both CBM-epu and CBM-npu result in one system failure in the $20 \times 5 = 100$ system runs,
whereas the two age-based policies each result in ten system failures.
System runtimes are less variable for the age-based policies
since they do not use the current system structure information, %available information on the component level,
i.e., which components have failed and which not. %during operational cycles
System runtimes are more variable for the condition-based policies,
but longer on average, leading to lower realized cost rates.
Again, CBM-cpu performs slightly better than CBM-epu
as it has the possibility to learn within an operational cycle as well,
managing to avoid system failures entirely. 
Still, CBM-npu does surprisingly well, considering that it never updates prior component parameters. 
This confirms the observation from Section~\ref{sec:ex-3}
that current system structure information delivers a major part of the gains for our CBM policy.

\begin{figure}
\includegraphics[width=\textwidth]{br1sim3fig5Tt01ac}
\caption{Simulation study with 20 repetitions, where each repetition consists of 5 subsequent operational cycles.
Case C: $\lambda_k = 2\ykz$.}
\label{fig:br1sim3fig5Tt01}
\end{figure}

\subsection{Case C: failure times later than expected, $\lambda_k = 2 \ykz$}
\label{sec:case3}

The third scenario consists of component failure times that are later than expected,
obtained by simulating Weibull lifetimes using $\lambda_k = 2 \ykz$ for all four component types.
The results are depicted in Figure~\ref{fig:br1sim3fig5Tt01}.

As the expert's information now underestimates component lifetimes,
all policies (except the corrective policy) act quite cautiously,
leading to complete avoidance of system failures.
However, the mean system runtimes $\mrsys$ exhibit
a clear hierarchy between age-based and condition-based policies,
and within these two groups a clear effect of the parameter update strategy.
Here, CBM-cpu obtains the longest system runtimes, and thus the lowest cost rates,
with a clear advantage over CBM-epu and CBM-npu.

Overall, it seems that CBM-cpu can exert its adaptive strengths most effectively
when prior assumptions are too pessimistic,
while also performing well when prior assumptions are correct or too optimistic.
%This means that experts may confidently tend towards pessimistic specifications
%of the component mean time to failures $\yzk$.


% ------------- old overview figures ----------------------------------
\iffalse
\begin{figure}
\includegraphics[width=0.5\textwidth]{br1sim1fig5}
\includegraphics[width=0.5\textwidth]{br1sim2fig5}
\includegraphics[width=0.5\textwidth]{br1sim3fig5}
\includegraphics[width=0.5\textwidth]{br1sim4fig5}
\caption{***Simulation study with ***prospective costrate criterion***, $\tthresh = 0.5$.
TL: $\lambda_k =     \ykz$,
TR: $\lambda_k = 0.5 \ykz$,
BL: $\lambda_k = 2   \ykz$,
BR: $\lambda_k = 0.2 \ykz$}
\label{fig:simprospcostrate05}
\end{figure}

\begin{figure}
\includegraphics[width=0.5\textwidth]{br1sim1fig5T}
\includegraphics[width=0.5\textwidth]{br1sim2fig5T}
\includegraphics[width=0.5\textwidth]{br1sim3fig5T}
\includegraphics[width=0.5\textwidth]{br1sim4fig5T}
\caption{***Simulation study with ***total costrate criterion***, $\tthresh = 0.5$.
TL: $\lambda_k =     \ykz$,
TR: $\lambda_k = 0.5 \ykz$,
BL: $\lambda_k = 2   \ykz$,
BR: $\lambda_k = 0.2 \ykz$}
\label{fig:simtotalcostrate05}
\end{figure}

\begin{figure}
\includegraphics[width=0.5\textwidth]{br1sim2fig5Tt01}
\includegraphics[width=0.5\textwidth]{br1sim2fig5Tt02}
\caption{***Simulation study with ***total costrate criterion***, early failures ($0.5 \times$ MTTF),
$\tthresh = 0.1$ (left) and $\tthresh = 0.2$ (right).}
\label{fig:simtotalcostrate0102}
\end{figure}
\fi


\section{Summary and outlook}
\label{sec:outlook}

In this paper, we have proposed a new condition-based maintenance policy
for systems based on monitoring of the functioning of components.
We use the survival signature to compute system reliability functions
for given component models, and propose to use Weibull component models
with conjugate priors to express expert knowledge.
The resulting policy shows promising results for
different numerical examples and a small simulation study.

There are several aspects in our modeling that warrant further investigation.
In this paper, we assume that component failure times are observed precisely via monitoring.
However, in many applications, the component status can only be determined
at inspections, leading to interval-censored failure time observations.
Extending the model to account for interval-censored failure times
would allow to derive a condition-based maintenance policy with wider-spaced inspections
by determining the optimal moment for the next inspection of the system.
%which could also determine the optimal moment for the next inspection of the system.
%*** extend the model for inspection-based CBM (if no repair necessary, find optimal timing of next inspection),
%we would need to add interval-censored data (failure has happened between last and current inspection)

Another way to relax the assumption of precise component monitoring
is to consider unreliable detection of component failures.
Our approach may be quite sensitive to the non-detection of component failures,
and one could tackle this by including a probabilistic observation model,
which assigns probabilities to false positives and false negatives in the component monitoring.

We also assume that the component-specific Weibull shape parameters are known,
and update only the scale parameters.
We restricted ourselves to this simplistic modelling
to allow for a clear and straight-forward demonstration of the general approach.
In a practical application, the shape parameters should be jointly updated together with the scale parameters.
This could be implemented by adopting the discretized approach by \cite{1969:soland}.

As mentioned at the end of Section~\ref{sec:costrate},
we currently disregard the possibility of common cause failures,
i.e., simultaneous failures with a shared root cause.
Common cause failures can drastically reduce the system reliability
by leading to sudden system failures.
To offer a realistic estimation of the optimal time to maintenance $\tausnow$,
the system RLD $\Rsys(t)$ needs to account for such failures when they are possible.
%This needs to be implemented on the level of component models,
%specifically $P\Big( \bigcap_{k=1}^K \{ C^k_t = l_k\} \Big)$ 
%(not a product anymore if common-case failures can affect several component types at once).
When common cause failures can affect several component types at once,
the joint posterior predictive probabilities $P\Big( \bigcap_{k=1}^K \{ C^k_t = l_k\} \Big)$ in \eqref{eq:sysrel-survsign}
cannot be obtained as a product of component posterior predictive probabilities any more.
We think that by combining ideas from \cite{2015:coolen-coolen-maturi} and \cite{Troffaes2014a},
the proposed policy could be extended to accommodate common cause failures.

In our policy, we assumed that in both preventive and corrective replacement
all components are restored to an as good as new state.
An obvious alternative is to only replace the failed components,
but would this offer any advantage?
Furthermore, could it be even optimal to replace a subset of the failed components only?
Such selective component replacement schemes
seem to be a promising area for research.
To account for different ages of components at system startup,
one could introduce artificial component types in \eqref{eq:sysrel-survsign},
grouping components both by type and by age.
However, computational complexity increases with the number of component types,
%and the conceptual algorithmic advantage of the survival signature vanishes
%when each component forms its own type
so it seems most promising to replace batches of components,
or all components of one type.
%*** study other replacement schemes, e.g., replacing only failed components:
%when repaired/replaced components are present in system,
%then one needs to use actual ages (and not time since system startup) in parameter update,
%introduces new component type for $\Rsysnow(t)$ calculation but with same posterior parameters as unreplaced component type.
%Repaired/replaced components require the creation of a separate type in the survival signature decomposition
%and thus the calculation of separate posterior predictive probabilities
%$P(C^k_t = l_k \mid \nkz,\ykz,\vectknow)$,
%which use, however, the same hyperparameter learning (using actual component ages) as the `parent type'.

A further simplifying assumption relates to the cost of preventive and corrective replacements,
which we consider as fixed.
In many real-life applications, maintenance costs vary over time
(e.g., increased labour cost due to overtime or night shift premiums),
and may depend on the number of failed components, or even on the exact system state.
In our approach, time-dependent cost structures can be easily accounted for,
by letting $c_p$ and $c_c$ in \eqref{eq:gtau} depend on $\tau$.
It will be more complex to account for costs that are dependent on the system state,
but our approach already provides the necessary information
in form of the posterior predictive distributions derived in Section~\ref{sec:postpred}.
%since the calculation of $\Rsysnow(t)$ involves  

Another promising extension would be to use sets of conjugate priors for the component models.
This approach is described in \cite{2016:walter-coolen},
leading to sets of residual life distributions $\Rsysnow(t)$.
Consequently, we would obtain a set of expected operational cycle cost rate functions $\gnow(\tau)$.
%and this would lead to sets of expected operational cycle cost rate functions $\gnow(\tau)$.
The determination of the (set of) optimal time(s) to maintenance $\tausnow$
is then non-trivial,
and would involve decision criteria discussed in the imprecise probability literature,
like E-admissibility, maximality, or maximin \citep[see, e.g., \S 8][]{itip}.


\section*{Acknowledgments}

Gero Walter was supported by the DINALOG project CAMPI
(``Coordinated Advanced Maintenance and Logistics Planning for the Process Industries'').
%
Both authors would like to thank Frank Coolen for inspiring discussions
and valuable hints on maintenance optimization criteria.

\section*{Bibliography}

\bibliographystyle{elsarticle-harv}

\bibliography{refs}

\end{document}