Robust Optimum Life-Testing Plans under Progressive Type-I Interval Censoring Schemes with Cost Constraint

Abstract

This paper considers optimal design problems for the Weibull distribution, which can be used to model symmetrical or asymmetrical data, in the presence of progressive interval censoring in life-testing experiments. Two robust approaches, Bayesian and minimax, are proposed to deal with the dependence of the D-optimality and c-optimality on the unknown model parameters. Meanwhile, the compound design method is applied to ensure a compromise between the precision of estimation of the model parameters and the precision of estimation of the quantiles. Furthermore, to make the design become more practical, the cost constraints are taken into account in constructing the optimal designs. Two algorithms are provided for finding the robust optimal solutions. A simulated example and a real life example are given to illustrate the proposed methods. The sensitivity analysis is also studied. These new design methods can help the engineers to obtain robust optimal designs for the censored life-testing experiments.

1. Introduction

Progressive censoring is frequently employed in life-testing experiments because it permits removing the test units at the points other than the final termination point, which can save experiment time and/or cost. Recently, the research on progressive censoring has grown very fast. For the relevant research progress one may refer to two important monographs (Balakrishnan et al. [1]; Balakrishnan and Cramer [2]) and the review article (Balakrishnan [3]). In applications, progressive type-I and type-II censoring are two important types of progressive censoring schemes. They are usually required to continuously observe the testing process under a given censoring scheme. However, due to the high cost and/or possible danger, it is sometimes infeasible to carry out continuous inspection in monitoring the test. Alternatively, an interval inspection scheme can be used, where only the number of failures between two consecutive inspections is recorded. Combining the concepts of the progressive censoring and the interval censoring, Aggarwala [4] developed a progressive type-I interval censoring (PIC-I) scheme.

Since Aggarwala [4], many scholars have studied the statistical inference for PIC-I data under various life distributions by the maximum likelihood method and/or Bayesian method. Some of them are Xiang and Tse [5], Ng and Wang [6], Chen and Lio [7], Lin and Lio [8], Singh and Tripathi [9], Lodhi and Tripathi [10], Budhiraja and Pradhan [11], Wu and Chang [12], and Alotaibi et al. [13]. In addition to statistical inference, many researchers have been focused on optimally obtaining PIC-I test plans. For more information on this direction one can refer to, e.g., Wu et al. [14], Lin et al. [15], Lin et al. [16], Singh and Tripathi [9], Roy and Pradhan [17], Budhiraja and Pradhan [18], Roy and Pradhan [19], and Wu et al. [20].

In most of the aforementioned studies on designing PIC-I test schemes, the model parameters that appear in the design criteria are assumed to be known, usually gained from previous studies. Designs obtained under these planning values are called locally optimal designs. If these values are uncertain, which is usually the case, or are incorrectly specified, then these designs may not be optimal. Therefore, it would be helpful for researchers to obtain the optimal designs, which are robust against misspecifications of the values of the model parameters. For the PIC-I life tests, Roy and Pradhan [17,19] proposed to employ proper prior distributions over the entire parameter space to describe the knowledge about the parameters. Then the designs obtained were called Bayesian optimal designs. Unfortunately, we sometimes do not have enough information to construct such prior distributions. In this case, the minimax method can be used to obtain robust designs. The resulting design is often referred to as the minimax optimal design. In fact, in the literature of experimental design, the Bayesian and minimax strategies have been widely used to overcome the dependency of the optimal design on the unknown parameters (see, e.g., Chaloner and Verdinelli [21], Atkinson, et al. [22], and Yue and Zhou [23]). However, in our knowledge, these problems have not been properly addressed when designing the reliability test.

Previous reviewed studies mainly focused on single design objective, such as estimation of model parameters (e.g., Wu et al. [14]) or the qth quantile of the life distribution or minimization of the total cost of life testing (e.g., Roy and Pradhan [19]). However, sometimes researchers prefer to achieve a design that meets multiple objectives at the same time. This design is called a multiple-objective design, which is robust with respect to the design objectives. In the accelerated life test (ALT), Pan and Yang [24] proposed a compound design criterion to obtain a dual-objective optimal design, which could make a trade-off between the model parameter estimation and the model-based prediction. Bhattacharya et al. [25] gave a new design method based on multi-criteria, taking variance and cost factors into account in the context of the hybrid censored life-testing experiment. Though there are some studies on robust designs in reliability life tests for multiple design objectives, robust optimum designs for PIC-I tests have received little attention in the literature so far.

In this paper, we study robust optimum designs for PIC-I life tests by using the Bayesian and minimax strategies to deal with the model parameters dependency problem and the compound design method to fulfill multiple design objectives. Two robust optimum design criteria and two algorithms are provided to calculate robust optimum designs when the experimental cost is taken into account. Our methods can help practitioners have easier access to robust optimum designs when the model parameters are uncertain and there are more than one design objective, and can be easily extended to obtain robust designs for other test schemes.

The rest of this paper is organized as follows. Section 2 introduces the PIC-I test plan and derives the Fisher information matrix (FIM). Section 3 provides definitions of the Bayesian compound and the minimax compound optimality criteria for the PIC-I test plans based on both D-optimality and c-optimality, termed $BcD$-optimality and $McD$-optimality, respectively. Section 4 presents two optimization algorithms, i.e., mixed-integer nonlinear optimization (MNO) and Particle swarm optimization (PSO), to derive optimal equal spaced PIC-I test plans and general PIC-I test plans, respectively. Numerical results are given in Section 5 and Section 6 to illustrate the proposed methods. Conclusions and discussions are made in Section 7.

2. Preliminaries

Suppose that the lifetime T of a test product follows the Weibull distribution with the probability density function (pdf)

$$\begin{array}{c}\hfill f(t;\eta ,\nu )={\displaystyle \frac{\nu }{{\eta }^{\nu }}}{t}^{\nu -1}exp\left\{-{\left({\displaystyle \frac{t}{\eta }}\right)}^{\nu }\right\},t>0,\eta ,\nu >0,\end{array}$$

(1)

and the cumulative distribution function (cdf)

$$\begin{array}{c}\hfill F(t;\eta ,\nu )=1-exp\left\{-{\left({\displaystyle \frac{t}{\eta }}\right)}^{\nu }\right\}.\end{array}$$

(2)

Let $Y=log\left(T\right)$. We can convert the Weibull distribution to the extreme value (Gumbel) distribution given by

$$f(y;\mu ,\sigma )=\frac{1}{\sigma }exp\left\{\frac{y-\mu }{\sigma }-exp\left(\frac{y-\mu }{\sigma }\right)\right\},-\infty <\mu <\infty ,\sigma >0,$$

(3)

where $\mu =log\eta$ is the location parameter and $\sigma =1/\nu$ is the scale parameter. The corresponding cdf can be written as

$$G(y;\mu ,\sigma )=1-exp\left\{-exp\left({\displaystyle \frac{y-\mu }{\sigma }}\right)\right\}.$$

(4)

Assume that a PIC-I scheme is employed, where all N units are simultaneously placed on a life test at the beginning of the experiment, and interval inspections are conducted at time points $t_1,t_2,\dots ,t_k$. At the jth inspection time $t_j$, $n_j$ failed units are observed and $r_j$ surviving units are randomly removed from the experiment. Let $q_j$ denote the probability that a unit fails in the jth time interval given that the failure has not occurred in an earlier time interval, i.e.,

$$q_j=P(y_{j-1}\le Y\le y_j\mid Y\ge y_j)={\displaystyle \frac{G\left(y_j\right)-G\left(y_{j-1}\right)}{1-G\left(y_{j-1}\right)}}=1-exp(h_j-h_{j-1}),$$

(5)

where $h_0=0$, $h_j=-exp\left(z_j\right),j=2,\dots ,k$ with $z_j=(y_j-\mu )/\sigma$ and $y_j=log\left(t_j\right)$.

Under the PIC-I test plan, the distribution of the number of failed units $n_j$ is binomial, i.e.,

$$n_j\mid n_{j-1},\dots ,n_1,r_{j-1},\dots ,r_1\sim \mathrm{binomial}(m_j,q_j),$$

(6)

where $m_1=N$ and $m_j=N-{\sum }_{s=1}^{j-1}{n}_s-{\sum }_{s=1}^{j-1}{r}_s,j=2,\dots ,k$ is the number of non-removed surviving units at the beginning of the jth inspection. Note that the number of removal units $r_j$ is a random variable due to the randomness of the variable $n_j$, and its value can be computed through the predetermined percentages of the remaining survived units $p_j$ (with $p_k=1$). That is, $r_j=(m_j-n_j)p_j$. All data collected from the PIC-I test plan are denoted by $\mathcal{D}=\{n_j,r_j,j=1,\dots ,k\}$.

Let $\mathbf{\theta }={(\mu ,\sigma )}^T$ be the vector of the model parameters. Based on the data $\mathcal{D}$ and the cdf given in Equation (4), the likelihood function of $\mathbf{\theta }$ can be written as

$$\begin{array}{ccc}\hfill L\left(\mathbf{\theta }\right)& =& \prod _{j=1}^k\left(\genfrac{}{}{0pt}{}{m_j}{n_j}\right){\left[1-exp(h_j-h_{j-1})\right]}^{n_j}{\left[exp(h_j-h_{j-1})\right]}^{m_j-n_j}.\hfill \end{array}$$

Then, the corresponding log-likelihood function is

$$\ell \left(\mathbf{\theta }\right)=\sum _{j=1}^k\left\{log\left(\genfrac{}{}{0pt}{}{m_j}{n_j}\right)+n_{j}log\left[1-exp(h_j-h_{j-1})\right]+\left(m_j-n_j\right)(h_j-h_{j-1})\right\}.$$

(7)

Therefore, the maximum likelihood estimates (MLEs) of the parameters $\mu$ and $\sigma$ can be obtained by solving the following likelihood equations:

$${\displaystyle \frac{\partial \ell \left(\mathbf{\theta }\right)}{\partial \mu }}={\displaystyle \frac{\partial \ell \left(\mathbf{\theta }\right)}{\partial \sigma }}=0.$$

Furthermore, the FIM of the parameters $\mathbf{\theta }$ can be written as

$$\mathcal{I}\left(\mathbf{\theta }\right)=\left(\begin{array}{cc}{\mathcal{I}}_{{\mu }^2}& {\mathcal{I}}_{\mu \sigma }\\ {\mathcal{I}}_{\mu \sigma }& {\mathcal{I}}_{{\sigma }^2}\end{array}\right),$$

(8)

where

$$\begin{array}{ccc}\hfill {\mathcal{I}}_{{\mu }^2}& =& \frac{1}{{\sigma }^2}\sum _{j=1}^{k}E\left(m_j\right){\displaystyle \frac{{(h_j-h_{j-1})}^2}{1-exp(h_j-h_{j-1})}}exp(h_j-h_{j-1}),\hfill \\ \hfill {\mathcal{I}}_{{\sigma }^2}& =& \frac{1}{{\sigma }^2}\sum _{j=1}^{k}E\left(m_j\right){\displaystyle \frac{{(z_j{h}_j-z_{j-1}{h}_{j-1})}^2}{1-exp(h_j-h_{j-1})}}exp(h_j-h_{j-1}),\hfill \\ \hfill {\mathcal{I}}_{\mu \sigma }& =& \frac{1}{{\sigma }^2}\sum _{j=1}^{k}E\left(m_j\right){\displaystyle \frac{(h_j-h_{j-1})(z_j{h}_j-z_{j-1}{h}_{j-1})}{1-exp(h_j-h_{j-1})}}exp(h_j-h_{j-1})\hfill \end{array}$$

and $E\left(m_j\right)=N{S}_{j-1},j=2,\dots ,k$. Here, $S_1=1$, and $S_{j-1}={\prod }_{s=1}^{j-1}(1-q_s)(1-p_s)$ is the survival probability of the test unit until the time point $t_{j-1}$. Under some mild regularity conditions, the FIM is approximate to the inverse of the variance-covariance matrix of the MLE of $\mathbf{\theta }$.

3. Robust PIC-I Test Plan

3.1. PIC-I Test Plan

In this section, we investigate the optimal PIC-I test plan with limited cost constraint. Generally, the  PIC-I test plan consists of the total number of test units N, the number of inspections k, the inspection time points $t_j,j=1,\dots ,k$ and the proportions of the removals at each time point $p_j,j=1,\dots ,k-1$. For convenience, we denote the PIC-I test plan as $\xi$, i.e.,

$$\xi =\left\{N,t_j,p_j,j=1,\dots ,k,p_k=1\right\}.$$

From the expression of the FIM given in (8), it can be easily concluded that if there is no limit on the removal proportions $p_j$ and if the other conditions are fixed, the optimal choices of $p_j(j=1,\dots ,k-1)$ are zeros. Some evidences can be found in the work by Roy and Pradhan [17,19]. However, this is not consistent with the topic discussed in this paper. To investigate the influence of the removal proportions $p_j(j=1,\dots ,k)$ on the optimal PIC-I test plan, we assume that the scheme of removals at each points is predetermined. Then, the PIC-I test plan can be re-expressed as

$$\xi =\left\{N,t_j;p_j,j=1,\dots ,k,p_k=1,N\in {\mathcal{N}}^{+},t_j\in {\mathcal{R}}^{+},p_j\in \mathcal{P}\right\},$$

(9)

where ${\mathcal{N}}^{+}$ is the set of positive integers, ${\mathcal{R}}^{+}$ is the set of positive natural number, $\mathcal{P}$ is the set of possible values of $p_j$ and will be predetermined by the experimenter.

For practical convenience, many studies assume that the inspection time intervals have equal lengths, say $\tau$. Then, the inspection time points can be re-expressed as $t_j=j\tau ,j=1,\dots ,k,k\in \mathcal{K},\mathcal{K}=\{1,2,\dots ,k^{\diamond }\}$ where $k^{\diamond }$ is the maximum number of inspections. In addition, the removal proportions at each inspection times can be assumed constant, i.e., $p_1=p_2=\dots =p_{k-1}=p\in \mathcal{P}$. Therefore, the design (9) can be reduced to

$$\xi =\left\{N,k,\tau ;p,N\in {\mathcal{N}}^{+},k\in \mathcal{K},\tau \in {\mathcal{R}}^{+},p\in \mathcal{P}\right\}.$$

(10)

Since the FIM $\mathcal{I}\left(\mathbf{\theta }\right)$ given in (8) depends not only on the parameter vector $\mathbf{\theta }$, but also the design $\xi$, we denote $\mathcal{I}\left(\mathbf{\theta }\right)$ by $\mathcal{I}(\xi ;\mathbf{\theta })$ in what follows.

3.2. D- and c-Optimal Design Criteria

The first design criterion we considered is the D-optimality for estimating the model parameter vector $\mathbf{\theta }$ as efficiently as possible. The  D-optimal criterion is defined as follows:

$${\mathsf{\Psi }}_D(\xi ;\mathbf{\theta })=-\frac{1}{2}log\mid \mathcal{I}(\xi ;\mathbf{\theta })\mid =-\frac{1}{2}log({\mathcal{I}}_{{\mu }^2}{\mathcal{I}}_{{\sigma }^2}-{\mathcal{I}}_{\mu \sigma }^2),$$

(11)

where ${\mathcal{I}}_{{\mu }^2},{\mathcal{I}}_{{\sigma }^2},{\mathcal{I}}_{\mu \sigma }$ are given in Equation (8).

Let $\mathsf{\Xi }$ be the set consisting of all possible designs $\xi$ in the form of (9) or (10). The design ${\xi }_D^{\ast }$, which satisfies ${\xi }_D^{\ast }={min}_{\xi \in \mathsf{\Xi }}{\mathsf{\Psi }}_D(\xi ;\mathbf{\theta })$ for given $\mathbf{\theta }$ is called a D-optimal design. The rationality of taking use of this design criterion can be found in Atkinson et al. [22] and has been discussed by many authors (e.g., Wu et al. [14], and Roy and Pradhan [17]).

In addition to the estimation of the model parameters $\mu$ and $\sigma$, an experimenter may also be interested in estimating the qth quantile lifetime of the units (Roy and Pradhan [19]). Let $y_q$ be the logarithm of the qth quantile lifetime of the units, and ${\mathbf{c}}_q={(1,c_q)}^T={(1,log(-log(1-q)))}^T$. Under the distribution (4), then $y_q$ can be written as $y_q={\mathbf{c}}_q^T\mathbf{\theta }$. By the invariance property of the MLE and the delta method, the distribution of the estimator ${\widehat{y}}_q$ is approximately normal with mean $y_q$ and variance ${\mathbf{c}}_q^T{\mathcal{I}}^{-1}(\xi ;\mathbf{\theta }){\mathbf{c}}_q$. To efficiently estimate $y_q$, the following c-optimality criterion is usually adopted:

$${\mathsf{\Psi }}_c(\xi ;\mathbf{\theta })=log\left[{\mathbf{c}}_q^T{\mathcal{I}}^{-1}(\xi ;\mathbf{\theta }){\mathbf{c}}_q\right]=log\left[({\mathcal{I}}_{{\mu }^2}-2c_q{\mathcal{I}}_{\mu \sigma }+c_q^2{\mathcal{I}}_{{\sigma }^2})/({\mathcal{I}}_{{\mu }^2}{\mathcal{I}}_{{\sigma }^2}-{\mathcal{I}}_{\mu \sigma }^2)\right].$$

(12)

A design ${\xi }_c^{\ast }$ which minimizes ${\mathsf{\Psi }}_c(\xi ;\mathbf{\theta })$ for given $\mathbf{\theta }$ and q over the design space $\mathsf{\Xi }$ is called a c-optimal design.

Remark 1.

The design criterion given in Equation (12) depends on the setting of q, which will be determined based on some practical consideration. To overcome the dependence, a nonnegative weight function $\mathrm{W}\left(q\right)$ satisfying ${\int }_0^1\mathrm{W}\left(q\right)dq=1$ can be used following the idea of Kundu [26]. Then, the design criterion will be

$${\mathsf{\Psi }}_c(\xi ;\mathbf{\theta })=log\left[{\int }_0^1({\mathcal{I}}_{{\mu }^2}-2c_q{\mathcal{I}}_{\mu \sigma }+c_q^2{\mathcal{I}}_{{\sigma }^2})/({\mathcal{I}}_{{\mu }^2}{\mathcal{I}}_{{\sigma }^2}-{\mathcal{I}}_{\mu \sigma }^2)\mathrm{W}\left(q\right)dq\right].$$

To compare the optimal PIC-I test plan ${\xi }_L^{\ast }$ with another arbitrary test plan $\xi$ under a given L-optimality criterion ($L\in \{D,c\}$), we define the following efficiency function

$${\mathrm{Eff}}_L({\xi }_L;·)=exp({\mathsf{\Psi }}_L({\xi }_L^{\ast };·)-{\mathsf{\Psi }}_L(\xi ;·)).$$

(13)

3.3. Bayesian Compound Design Criterion

Note that the design criteria ${\mathsf{\Psi }}_D(\xi ;\mathbf{\theta })$ and ${\mathsf{\Psi }}_c(\xi ;\mathbf{\theta })$ depend not only on the design $\xi$, but also on the the parameters $\mathbf{\theta }$. Optimal designs obtained under a perfect guess (planning) values ${\mathbf{\theta }}_t$ are called locally optimal designs (Chernoff [27]). Numerical results given in Wu et al. [14] indicate that D-optimal test plans depend on the setting of the parameters $\mu$ and $\sigma$. To reduce the risks caused by misspecifying the planning values ${\mathbf{\theta }}_t$, many approaches, such as the Bayesian, minimax, adaptive, or sequential approaches have been proposed in the literature of optimal experimental designs, see Atkinson et al. [22] for details. However, planning robust test schemes has rarely been found in a reliability study. Therefore, in the following part of this subsection, we apply robust design techniques to obtain robust PIC-I test plans. To be consistent with most experiments in practical applications, we confine ourselves on the static robust design methods, including the Bayesian and the minimax approaches to obtain robust test plans against the uncertainty of the parameters $\mathbf{\theta }$.

In addition, even when all parameters in the D-optimality and c-optimality criteria have the same settings, the D-optimal design ${\xi }_D^{\ast }$ and the c-optimal design ${\xi }_c^{\ast }$ are not necessarily identical, see the numerical results listed in

Table 1. Locally optimal PIC-I test plans for different removal proportions p, planning values of $\mathbf{\theta }$, and design criteria.

p	$\mathit{\mu }$	$\mathit{\sigma }$	D-Optimality						c-Optimality
			$\mathit{N}$	$\mathit{k}$	$\tau$	$\mathit{k}\tau$	${\mathsf{\Psi }}_{\mathit{D}}$	${\mathrm{Eff}}_{\mathit{D}}\left({\xi }_{\mathit{c}}^{\ast }\right)$	$\mathit{N}$	$\mathit{k}$	$\tau$	$\mathit{k}\tau$	${\mathsf{\Psi }}_{\mathit{c}}$	${\mathrm{Eff}}_{\mathit{c}}\left({\xi }_{\mathit{D}}^{\ast }\right)$
0.1	$\widehat{\mu }-s_{\widehat{\mu }}$	$\widehat{\sigma }-s_{\widehat{\sigma }}$	74	6	2.0121	12.0724	−6.1284	0.9079	74	4	1.9860	7.9439	−4.0028	0.9975
0.1	$\widehat{\mu }-s_{\widehat{\mu }}$	$\widehat{\sigma }+s_{\widehat{\sigma }}$	74	9	2.1483	19.3346	−5.1970	0.9272	74	6	1.7824	10.6944	−3.1115	0.9882
0.1	$\widehat{\mu }+s_{\widehat{\mu }}$	$\widehat{\sigma }-s_{\widehat{\sigma }}$	74	6	2.6201	15.7203	−6.1268	0.9076	74	4	2.5874	10.3497	−4.0018	0.9970
0.1	$\widehat{\mu }+s_{\widehat{\mu }}$	$\widehat{\sigma }+s_{\widehat{\sigma }}$	73	9	2.7913	25.1220	−5.1945	0.8753	74	5	2.4145	12.0724	−3.1102	0.9873
0.3	$\widehat{\mu }-s_{\widehat{\mu }}$	$\widehat{\sigma }-s_{\widehat{\sigma }}$	74	4	3.0280	12.1120	−5.8608	0.9886	74	3	3.0755	9.2264	−3.8758	0.9981
0.3	$\widehat{\mu }-s_{\widehat{\mu }}$	$\widehat{\sigma }+s_{\widehat{\sigma }}$	74	6	3.1140	18.6839	−4.9051	0.9493	74	4	2.6947	10.7789	−3.0115	0.9835
0.3	$\widehat{\mu }+s_{\widehat{\mu }}$	$\widehat{\sigma }-s_{\widehat{\sigma }}$	74	4	3.9453	15.7813	−5.8592	0.9888	74	3	4.0074	12.0222	−3.8746	0.9977
0.3	$\widehat{\mu }+s_{\widehat{\mu }}$	$\widehat{\sigma }+s_{\widehat{\sigma }}$	74	6	4.0537	24.3223	−4.9027	0.9499	74	4	3.5113	14.0454	−3.0101	0.9828
0.1	$\widehat{\mu }-2s_{\widehat{\mu }}$	$\widehat{\sigma }-2s_{\widehat{\sigma }}$	74	6	1.4656	8.7933	−6.7534	0.9369	74	4	1.5798	6.3191	−4.6350	0.9969
0.1	$\widehat{\mu }-2s_{\widehat{\mu }}$	$\widehat{\sigma }+2s_{\widehat{\sigma }}$	73	11	1.8835	20.7189	−4.8260	0.8784	74	6	1.4889	8.9334	−2.7718	0.9797
0.1	$\widehat{\mu }+2s_{\widehat{\mu }}$	$\widehat{\sigma }-2s_{\widehat{\sigma }}$	74	6	2.4859	14.9152	−6.7508	0.9361	74	4	2.6829	10.7315	−4.6332	0.9961
0.1	$\widehat{\mu }+2s_{\widehat{\mu }}$	$\widehat{\sigma }+2s_{\widehat{\sigma }}$	73	10	3.2318	32.3185	−4.8203	0.8794	74	6	2.5223	15.1340	−2.7691	0.9770
0.3	$\widehat{\mu }-2s_{\widehat{\mu }}$	$\widehat{\sigma }-2s_{\widehat{\sigma }}$	74	3	3.1648	9.4943	−6.4923	0.3398	74	2	4.2754	8.5508	−4.4716	0.9275
0.3	$\widehat{\mu }-2s_{\widehat{\mu }}$	$\widehat{\sigma }+2s_{\widehat{\sigma }}$	74	7	2.7790	19.4529	−4.5149	0.9059	74	4	2.2228	8.8912	−2.6749	0.9630
0.3	$\widehat{\mu }+2s_{\widehat{\mu }}$	$\widehat{\sigma }-2s_{\widehat{\sigma }}$	74	2	5.4124	10.8248	−6.4900	0.3398	74	2	7.2750	14.5501	−4.4691	0.9287
0.3	$\widehat{\mu }+2s_{\widehat{\mu }}$	$\widehat{\sigma }+2s_{\widehat{\sigma }}$	73	6	4.7729	28.6373	−4.5095	0.9078	74	4	3.7753	15.1014	−2.6723	0.9595

. To get a good test plan that satisfies many design objectives, the compound design criterion (Cook and Wong [28]) should be a good choice. It can ensure that the optimal design has good performances for different design purposes (Pan and Yang [24]). Therefore, considering the uncertainty of the parameters in the model and two possible design objectives, we propose the Bayesian compound optimality criterion (termed $BcD$-optimality), which is given as follows:

$${\mathsf{\Psi }}_{BcD}(\xi ;\kappa )=\kappa {\int }_{\mathsf{\Theta }}{\mathsf{\Psi }}_D(\xi ;\mathbf{\theta })p\left(\mathbf{\theta }\right)d\mathbf{\theta }+(1-\kappa ){\int }_{\mathsf{\Theta }}{\mathsf{\Psi }}_c(\xi ;\mathbf{\theta })p\left(\mathbf{\theta }\right)d\mathbf{\theta },$$

(14)

where $p\left(\mathbf{\theta }\right)$ is the prior over the set $\mathsf{\Theta }$, which will be specified later, and $\kappa$ is the weight parameter reflecting the relative importance between the D-optimality and the c-optimality. A design, ${\xi }_{BcD}^{\ast }$, which minimizes ${\mathsf{\Psi }}_{BcD}(\xi ;\kappa )$ over the design space $\mathsf{\Xi }$ is called a $BcD$-optimal design. In this paper, we assume that the prior of $\mu$ is a censored normal distribution over the interval $[a,b]$, and its pdf is

$$\begin{array}{ccc}\hfill \pi (\mu ;{\mu }_0,{\sigma }_0^2)& =& {\displaystyle \frac{f(\mu ;{\mu }_0,{\sigma }_0^2)}{\mathsf{\Phi }(b;{\mu }_0,{\sigma }_0^2)-\mathsf{\Phi }(a;{\mu }_0,{\sigma }_0^2)}},a\le \mu \le b,{\sigma }_0>0,\hfill \end{array}$$

(15)

where $f(\mu ;{\mu }_0,{\sigma }_0^2)$ and $\mathsf{\Phi }(·;{\mu }_0,{\sigma }_0)$ are the pdf and cdf of the normal distribution $N({\mu }_0,{\sigma }_0^2)$ and ${\mu }_0,{\sigma }_0^2$ are the hyperparameters predetermined by the experimenter. The prior of the scale parameter $\sigma$ is assumed to be a censored inverse $\mathsf{\Gamma }({\nu }_0,{\gamma }_0)$ over the interval $[c,d]$ with the pdf

$$\begin{array}{ccc}\hfill \pi (\sigma ;{\nu }_0,{\gamma }_0)& =& {\displaystyle \frac{f(\sigma ;{\nu }_0,{\gamma }_0)}{{\mathsf{\Gamma }}^{-1}(d;{\nu }_0,{\gamma }_0)-{\mathsf{\Gamma }}^{-1}(c;{\nu }_0,{\gamma }_0)}},\hfill \\ \hfill f(\sigma ;{\nu }_0,{\gamma }_0)& =& {\displaystyle \frac{\mathsf{\Gamma }\left({\nu }_0\right)}{{\gamma }_0^{{\nu }_0}}}{\sigma }^{-({\nu }_0+1)}exp\left(-\frac{{\gamma }_0}{\sigma }\right),{\nu }_0,{\gamma }_0>0,\hfill \end{array}$$

where $f(·;{\nu }_0,{\gamma }_0)$ and ${\mathsf{\Gamma }}^{-1}(·;{\nu }_0,{\gamma }_0)$ are the pdf and cdf of the inverse Gamma distribution ${\mathsf{\Gamma }}^{-1}({\nu }_0,{\gamma }_0)$, respectively, and  ${\nu }_0,{\gamma }_0$ are the hyperparameters predetermined by the experimenter. Under the assumption of the priors of $\mu$ and $\sigma$ being independent, the joint prior density of $\mu$ and $\sigma$ is given as

$$\pi (\mu ,\sigma )=\pi (\mu ;{\mu }_0,{\sigma }_0^2)·\pi (\sigma ;{\nu }_0,{\gamma }_0),(\mu ,\sigma )\in \mathsf{\Theta }=[a,b]\times [c,d].$$

In addition, when the weight $\kappa$ equals 0 (or 1), then the $BcD$-optimal design ${\xi }_{BcD}^{\ast }$ will reduce to Bayesian c (or D)-optimal design and be denoted as ${\xi }_{Bc}^{\ast }$ (or ${\xi }_{BD}^{\ast }$). In the special case that the prior has only one support point ${\mathbf{\theta }}_t$, the $BcD$-optimal design becomes the locally $cD$-optimal design, which can be denoted as ${\xi }_{cD}^{\ast }$ and the corresponding design criterion is

$${\mathsf{\Psi }}_{cD}(\xi ;\kappa ,{\mathbf{\theta }}_t)=\kappa {\mathsf{\Psi }}_D(\xi ;{\mathbf{\theta }}_t)+(1-\kappa ){\mathsf{\Psi }}_c(\xi ;{\mathbf{\theta }}_t).$$

(16)

To deal with the integration in the Bayesian optimality criterion, we can draw Q samples ${\mathbf{\theta }}_l$ from the joint prior distribution $\pi (\mu ,\sigma )$ and use the MCMC method to obtain the approximation of the Bayesian optimality criterion (14) (Roy and Pradhan [17,19]), i.e.,

$${\mathsf{\Psi }}_{BcD}(\xi ;\kappa )=\frac{1}{Q}\sum _{l=1}^Q\left\{-\frac{\kappa }{2}log\mid \mathcal{I}(\xi ;{\mathbf{\theta }}_l)\mid +(1-\kappa )log\left[{\mathbf{c}}_q^T{\mathcal{I}}^{-1}(\xi ;{\mathbf{\theta }}_l){\mathbf{c}}_q\right]\right\}.$$

However, the computational burden for this approximation will increase rapidly with the sample size Q. To reduce computational burden, Foo and Duffull [29] proposed a hypercube D-optimality criterion to optimize the logarithm of the product of the normalized determinants of FIM over the set of their model parameters, which consists of the combinations of the 2.5th and 97.5th percentiles of the priors distributions. In this paper, we utilize a similar idea to alleviate the computational burden, but more percentiles of the prior parameter distribution are used to ensure the accuracy of approximation. We denote the set of all combinations of these percentiles by ${\mathsf{\Theta }}_Q$. Then

$${\mathsf{\Psi }}_{BcD}(\xi ;\kappa )=\frac{1}{Q}\sum _{{\mathbf{\theta }}_l\in {\mathsf{\Theta }}_Q}\left\{-\frac{\kappa }{2}log\mid \mathcal{I}(\xi ;{\mathbf{\theta }}_l)\mid +(1-\kappa )log\left[{\mathbf{c}}_q^T{\mathcal{I}}^{-1}(\xi ;{\mathbf{\theta }}_l){\mathbf{c}}_q\right]\right\}.$$

(17)

3.4. Minimax Compound PIC-I Test Plan

Another popular robust design approach is the minimax method, which aims at finding the design that minimizes the maximum values of the compound optimality criterion over the parameters space ${\mathsf{\Theta }}_Q$. Similar to the Bayesian optimality criterion, the minimax optimality criterion (termed $McD$-optimality) is defined as follows:

$${\mathsf{\Psi }}_{McD}(\xi ;\omega )=\omega \underset{\mathbf{\theta }\in {\mathsf{\Theta }}_Q}{max}{\mathsf{\Psi }}_D(\xi ;\mathbf{\theta })+(1-\omega )\underset{\mathbf{\theta }\in {\mathsf{\Theta }}_Q}{max}{\mathsf{\Psi }}_c(\xi ;\mathbf{\theta }),$$

(18)

where $\omega$ is the weight parameter similar as $\kappa$. The design ${\xi }_{McD}^{\ast }$ minimizing ${\mathsf{\Psi }}_{McD}(\xi ;\omega )$ over the design space $\mathsf{\Xi }$ is called a $McD$-optimal design. In addition, when the weight $\omega$ equals 0 (or 1), then the $McD$-optimal design ${\xi }_{McD}^{\ast }$ will reduce to the minimax c-optimal design and be denoted as ${\xi }_{Mc}^{\ast }$ (or ${\xi }_{MD}^{\ast }$).

3.5. Cost Constraint

To make the test plan be more practical for experimenters, we take the budget of the experiment into account at the planning stage. For the general design $\xi$ defined in Equation (9), we assume that there is a restriction on the total cost

$$C_1\left(\xi \right)=N{C}_s+k{C}_i+t_k{C}_o\le C,$$

(19)

where $C_s$ is the cost of one test unit, $C_i$ is the cost of one inspection, $C_o$ is the operation cost of one unit time, and C is the total cost. For more details one can refer to Wu et al. [14]. Furthermore, for the special test plan $\xi$ defined in (10), the cost constraint is defined as follows:

$$C_2\left(\xi \right)=N{C}_s+k{C}_i+k\tau C_o\le C.$$

(20)

4. Algorithms

4.1. Mixed-Integer Nonlinear Optimization Algorithm

In this subsection, we give an algorithm inspired by Wu et al. [14] to search for the robust PIC-I test plan with equal inspection intervals of length $\tau$ and the constant removal proportion p, except at the end of the experiment. The design is given in (10) and the cost constraint is presented in (20). For clarity, we rewrite the robust design criteria ${\mathsf{\Psi }}_{BcD}(\xi ;\kappa )$ and ${\mathsf{\Psi }}_{McD}(\xi ;\omega )$ by ${\mathsf{\Psi }}_{BcD}(N,k,\tau ;\kappa )$ and ${\mathsf{\Psi }}_{McD}(N,k,\tau ;\omega )$, respectively. Then, the optimization problem can be expressed as

$$\begin{array}{ccc}\hfill & & \min {\mathsf{\Psi }}_{BcD}(N,k,\tau ;\kappa )\mathrm{or}{\mathsf{\Psi }}_{McD}(N,k,\tau ;\omega )\hfill \\ & & \mathrm{s}.\mathrm{t}.C_2\left(\xi \right)=N{C}_s+k{C}_i+k\tau C_o\le C,\hfill \\ & & N\in {\mathcal{N}}^{+},\tau \in {\mathcal{R}}^{+},k\in \mathcal{K}.\hfill \end{array}$$

(21)

Since $C_2\left(\xi \right)\le C$, we have $N\le (C-(k{C}_i+k\tau C_o))/C_s\le \left(C-C_i\right)/C_s$. Thus, the upper bound of N is obtained. Because the D-optimality and c-optimality are decreasing functions of N, we substitute N by its upper bound $(C-(k{C}_i+k\tau C_o))/C_s$ in the design criteria. Furthermore, for a given value of N, we have $k\le (C-N{C}_s)/C_i\le (C-C_s)/C_i$. Then the upper bound of k is $min\{(C-C_s)/C_i,k^{\diamond }\}$. Similarly, we can obtain the upper bound of $\tau$, $\tau \le (C-C_s-k{C}_i)/\left(k{C}_o\right)$ for a given k. Therefore, the algorithm to solve the optimization problem (21) is given in Algorithm 1.

Algorithm 1 MNO Algorithm.

1: Set the values of the cost parameters $C_s,C_i,C_o,C$; the maximum number of inspection $k^{\diamond }$; the removal proportion p; the region of the model related parameters $\mathsf{\Theta }$; the hyperparamters in the priors for Bayesian optimality criterion; the weight parameter $\kappa$ or $\omega$. 2: Compute the upper bound of the number of inspections $$\tilde{k}=min\left\{\left[\frac{C-C_s}{C_i}\right],k^{\diamond }\right\},$$ where $\left[y\right]$ denotes the greatest integer less than or equal to y. 3: Set $k=2$. 4: Calculate the upper bound of the length of inspection interval $${\tau }_k=\frac{C-C_s-k{C}_i}{k{C}_o}.$$ 5: For a given k, by  the optimization method, such as the PSO algorithm or the grid method, if the minimum length of the inspection interval is to be considered due to practical constraints, find the solution ${\tau }_k^{\ast }$ to the problem $${\tau }_k^{\ast }=arg\underset{\tau \in [0,{\tau }_k]}{min}{\mathsf{\Psi }}_{BcD}(N(k,\tau ),k,\tau ;\kappa )\mathrm{or}{\mathsf{\Psi }}_{McD}(N(k,\tau ),k,\tau ;\omega ),$$ where $N(k,\tau )=(C-(k{C}_i+k\tau C_o))/C_s$. 6: Set $k=k+1$. If $k\le \tilde{k}$ go to Step 4, else go to Step 7. 7: Find the optimal design ${\xi }^{\ast }=\left\{N^{\ast },k^{\ast },{\tau }^{\ast }\right\}$ such that $$\begin{gathered} \begin{array}{c}{\mathsf{\Psi }}_{BcD}(N^{\ast },k^{\ast },{\tau }^{\ast };\kappa )=\underset{1\le k\le \tilde{k}}{min}{\mathsf{\Psi }}_{BcD}(N(k,{\tau }_k^{\ast }),k,{\tau }_k^{\ast };\kappa )\mathrm{or} \\ {\mathsf{\Psi }}_{McD}\left({\xi }^{\ast }\right)=\underset{1\le k\le \tilde{k}}{min}{\mathsf{\Psi }}_{McD}(N(k,{\tau }_k^{\ast }),k,{\tau }_k^{\ast };\omega ).\hfill \end{array} \end{gathered}$$ 8: Output the optimal test plan ${\xi }^{\ast }=\{N^{\ast },k^{\ast },{\tau }^{\ast }\}$.

4.2. Particle Swarm Optimization Algorithm

In Section 4.1, we have considered a special case of the PIC-I test plan and give an algorithm to obtain the robust test plan. However, the procedure depends on the cost function and the assumption of equal length of inspection interval and the same removal proportion. In this subsection, we give an algorithm to solve the general design problem defined as follows:

$$\begin{array}{ccc}\hfill & & \min {\mathsf{\Psi }}_{BcD}(\xi ;\kappa )\mathrm{or}{\mathsf{\Psi }}_{McD}(\xi ;\omega )\hfill \\ & & \mathrm{s}.\mathrm{t}.C_1\left(\xi \right)\le C,\hfill \\ & & N\in {\mathcal{N}}^{+},p_i\in \mathcal{P},0\le t_1\le t_2\le \dots \le t_k,t_i\in (0,t_{max}).\hfill \end{array}$$

(22)

The concrete form of the design $\xi$ has been given in (9). From the definition of the cost function $C_1\left(\xi \right)$ in (19), we have $N\le (C-k{C}_i-t_k{C}_o)/C_s$. The information given in (8) implies that more test units will provide more information. Therefore, we substitute N in (8) by the supreme value $(C-k{C}_i-t_k{C}_o)/C_s$ to solve the optimization problem (22). Then, the optimality criterion ${\mathsf{\Psi }}_{BcD}\left(\xi \right)$ or ${\mathsf{\Psi }}_{McD}\left(\xi \right)$ is a function of $p_i,t_i$, which are all continuous variables. To solve this continuous optimization problem with constraints, we use the PSO algorithm, which is a population based stochastic optimization method inspired by the social behavior of birds flocking or fish schooling. It was introduced by Eberhart and Kennedy [30], improved by many authors to deal with all kinds of optimization problems. It has shown high efficiency in finding optimal points in various disciplines (Poli et al. [31], Ruidas et al. [32,33]). The PSO algorithm is derivative-free. There are few tuning parameters required of the algorithm and the knowledge of good solutions is retained by all particles, and particles in the swarm share information between, which makes the algorithm easily escape from local minima and converge at a fast rate. Recently, different versions of PSO have been used to solve all kinds of optimal design problems (see Chen et al. [34], Zhou et al. [35], and Liu et al. [36]). In the following we give a summary of the PSO algorithm for completeness. For the sake of brevity and clarity, we vectorize the design $\xi$ and denote it by $X=[t_1,t_2,\dots ,t_k]$, and by $f\left(X\right)$ denote the optimality criterion function ${\mathsf{\Psi }}_{BcD}(\xi ;\kappa )$ or ${\mathsf{\Psi }}_{McD}(\xi ,\omega )$. The optimization problem (22) is rewritten as

$$\begin{array}{ccc}\hfill & & \min f\left(X\right)\hfill \\ & & \mathrm{s}.\mathrm{t}.C_2\left(X\right)-C\le 0,\hfill \\ & & 0\le t_1\le t_2\le \dots \le t_k,t_i\in (0,t_{max}),0<k<k^{\diamond }.\hfill \end{array}$$

(23)

Here, $u_i,i=1,2$ are real random numbers between 0 and 1, $P_i$ is the best candidate solution found for the ith particle, $P_g$ is the best candidate solution for the entire population particles, and $\alpha$, ${\alpha }_i,i=1,2$ are the user defined coefficients which respectively control the inertia, the exploitive, and the explorative attributes of the particle motion.

5. Numerical Example

In this section, we present a numerical example to illustrate the applications of the proposed robust design methods. As in the study by Wu et al. [14], we use the algorithm given by Aggarwala [4] to generate data with $n=20, k=5, \tau =2, \eta =5, \nu =2$ and the predetermined removed proportions $(p_1,\dots ,p_5)=(0,0.2,0.3,0.4,1)$. Then the corresponding parameters in model (4) become $\mu =log5 ,\sigma =1/2$. The generated data are presented in

Table 2. Progressively type-I interval-censored samples.

i	1	2	3	4	5
$n_i$	2	4	6	2	1
$r_i$	0	2	1	1	1
$p_i$	0	0.2	0.3	0.4	1

. Based on the data, we easily obtain the MLEs of $\mu$ and $\sigma$, $\widehat{\mu }=1.8454$ and $\widehat{\sigma }=0.5091$, and their standard errors $s_{\widehat{\mu }}=0.1329$ and $s_{\widehat{\sigma }}=0.1157$, respectively. To use the Bayesian or minimiax design criteria, we assume that the range of the values of $\mu$ and $\sigma$ are $[a,b]=[\widehat{\mu }-s_{\widehat{\mu }},\widehat{\mu }+s_{\widehat{\mu }}]$ and $[c,d]=[\widehat{\sigma }-s_{\widehat{\sigma }},\widehat{\sigma }+s_{\widehat{\sigma }}]$, respectively. The hyperparameters in the censored normal distribution $N({\mu }_0,{\sigma }_0^2)$ are ${\mu }_0=1.8$ and ${\sigma }_0=0.2$, respectively. Similarly, the hyperparameters in the censored inverse distribution ${\mathsf{\Gamma }}^{-1}({\nu }_0,{\gamma }_0)$ are assumed to be ${\nu }_0=27$ and ${\gamma }_0=13$, respectively, which implies that the mean of ${\mathsf{\Gamma }}^{-1}({\nu }_0,{\gamma }_0)$ is $0.5$ and the variance is $0.01$. The cost parameters are assumed as follows: $C_s=\$80$/unit, $C_i=\$3$, $C_o=\$2.5$/10 h, and $C=\$6000$. In our numerical results, the set ${\mathsf{\Theta }}_Q=\{({\mu }_{i\mathrm{th}},{\sigma }_{j\mathrm{th}}),i,j=0,1,\dots ,10\}$, where ${\mu }_{i\mathrm{th}}$ is the ith quantile of $\pi (\mu ;{\mu }_0,{\sigma }_0^2)$ and ${\sigma }_{j\mathrm{th}}$ is the jth quantile of $\pi (\sigma ;{\nu }_0,{\gamma }_0)$, which are used to obtain the approximation of the $BcD$-optimality criterion in (14) or to calculate the $McD$-optimality criterion in (18). We set $q=0.1$ in the $BcD$- and $McD$-optimal design criteria. Since the $McD$-optimal PIC-I test plans are very similar to their corresponding $BcD$-optimal ones, then we do not show the numerical $McD$-optimal PIC-I test plans in what follows to save space. We will give some concluding remarks in Section 6. The R codes, written to obtain the results in this paper, can be obtained from the authors upon request.

5.1. Locally Optimal PIC-I Test Plans

By Algorithm 1, we first compute the locally D- and c-optimal equal spaced (ES) PIC-I test plans (denoted as ${\xi }_D^{\ast }$ and ${\xi }_c^{\ast }$, respectively) at the planning values $\mu =log5$ and $\sigma =1/2$ when the removal proportions are $p=0.1$ and $p=0.3$, respectively.

$$\begin{array}{ccc}\hfill p=0.1:{\xi }_D^{\ast }:& & N=74,k=7,\tau =1.9261,{\mathsf{\Psi }}_D=-5.6620,\hfill \\ & & {\mathrm{Eff}}_D\left({\xi }_{BD}^{\ast }\right)=0.9834,{\mathrm{Eff}}_D\left({\xi }_{MD}^{\ast }\right)=0.9735,\hfill \\ \hfill {\xi }_c^{\ast }:& & N=74,k=5,\tau =1.7235,{\mathsf{\Psi }}_c=-3.5486,\hfill \\ & & {\mathrm{Eff}}_c\left({\xi }_{Bc}^{\ast }\right)=0.9872,{\mathrm{Eff}}_c\left({\xi }_{Mc}^{\ast }\right)=0.9911.\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill p=0.3:{\xi }_D^{\ast }:& & N=74,k=5,\tau =2.7647,{\mathsf{\Psi }}_D=-5.3891,\hfill \\ & & {\mathrm{Eff}}_D\left({\xi }_{BD}^{\ast }\right)=0.9609,{\mathrm{Eff}}_D\left({\xi }_{MD}^{\ast }\right)=0.9485,\hfill \\ & & {\xi }_c^{\ast }:N=74,k=3,\tau =2.6525,{\mathsf{\Psi }}_c=-3.4414,\hfill \\ & & {\mathrm{Eff}}_c\left({\xi }_{Bc}^{\ast }\right)=0.9622,{\mathrm{Eff}}_c\left({\xi }_{Mc}^{\ast }\right)=0.9824.\hfill \end{array}$$

(25)

In Equations (24) and (25), we also present the efficiencies of the Bayesian and minimax D- and c-optimal test plans with $\mathsf{\Theta }={\mathsf{\Theta }}_1$, which will be defined later. These efficiencies indicate that these robust designs perform well.

To assess the influence of the settings of the parameters $\mu$ and $\sigma$ and the design criteria, the optimal PIC-I test plans under different parameters settings and removal proportions for different design criteria are computed, which are shown in Table 1. Comparing with the optimal test plans presented in Equations (24) and (25), we find that the settings of the planning values will have a great impact on the final designs, especially when the uncertainty of $\mu$ and $\sigma$ encountered in the design stage becomes great—so do the removal proportions and design criteria. The efficiencies of the c- and D-optimal designs relative to the D- and c-optimal designs are also shown in Table 1. These efficiencies imply that the design purposes will have an effect on the finial designs. In addition, for some planning values, the c-optimal test plans can be very inefficient, whereas the D-optimal test plans perform well under the c-optimality in all cases considered here.

5.2. Bayesian Optimal PIC-I Test Plans

To obtain optimal PIC-I test plans that fulfill multiple design purposes, we compute $cD$-optimal compound designs for different weights $\kappa$, planning values $\mathbf{\theta }$, and removal proportions p. Some of the results are given in

Table 3. $cD$-optimal PIC-I test plans for different weights $\kappa$, different removal proportions p, and the different planning values $\mathbf{\theta }$.

p	$\mathit{\kappa }$	$\mathit{\mu }=log5,\mathit{\sigma }=0.5$							$\mathit{\mu }=\widehat{\mathit{\mu }}-2{\mathit{s}}_{\widehat{\mathit{\mu }}},\mathit{\sigma }=\widehat{\mathit{\sigma }}-2{\mathit{s}}_{\widehat{\mathit{\sigma }}}$
		$\mathit{N}$	$\mathit{k}$	$\tau$	$\mathit{k}\tau$	${\mathsf{\Psi }}_{\mathit{cD}}$	${\mathrm{Eff}}_{\mathit{D}}\left({\xi }_{\mathit{cD}}^{\ast }\right)$	${\mathrm{Eff}}_{\mathit{c}}\left({\xi }_{\mathit{cD}}^{\ast }\right)$	$\mathit{N}$	$\mathit{k}$	$\tau$	$\mathit{k}\tau$	${\mathsf{\Psi }}_{\mathit{cD}}$	${\mathrm{Eff}}_{\mathit{D}}\left({\xi }_{\mathit{cD}}^{\ast }\right)$	${\mathrm{Eff}}_{\mathit{c}}\left({\xi }_{\mathit{cD}}^{\ast }\right)$
0.1	0.0	74	5	1.7235	8.6174	−3.5486	0.9359	1.0000	74	4	1.5798	6.3191	−4.6350	0.9369	1.0000
	0.1	74	7	1.7453	12.2173	−3.7574	0.9961	0.9976	74	5	1.5582	7.7908	−4.8458	0.9981	0.9990
	0.2	74	7	1.7762	12.4334	−3.9686	0.9974	0.9973	74	5	1.5587	7.7933	−5.0575	0.9981	0.9990
	0.3	74	7	1.8025	12.6172	−4.1800	0.9983	0.9971	74	5	1.5590	7.7951	−5.2693	0.9981	0.9990
	0.4	74	7	1.8254	12.7781	−4.3915	0.9989	0.9967	74	5	1.5593	7.7964	−5.4810	0.9981	0.9990
	0.5	74	7	1.8460	12.9219	−4.6031	0.9993	0.9964	74	6	1.5021	9.0127	−5.6929	0.9998	0.9975
	0.6	74	7	1.8646	13.0523	−4.8148	0.9996	0.9960	74	6	1.4936	8.9619	−5.9049	0.9999	0.9974
	0.7	74	7	1.8817	13.1720	−5.0265	0.9998	0.9957	74	6	1.4858	8.9146	−6.1170	0.9999	0.9973
	0.8	74	7	1.8975	13.2828	−5.2383	0.9999	0.9953	74	6	1.4785	8.8709	−6.3291	1.0000	0.9972
	0.9	74	7	1.9123	13.3858	−5.4501	1.0000	0.9949	74	6	1.4717	8.8305	−6.5413	1.0000	0.9970
	1.0	74	7	1.9261	13.4827	−5.6620	1.0000	0.9946	74	6	1.4656	8.7935	−6.7534	1.0000	0.9969
0.3	0.0	74	3	2.6524	7.9573	−3.4414	0.9312	1.0000	74	2	4.2754	8.5508	−4.4716	0.3398	1.0000
	0.1	74	4	2.6559	10.6235	−3.6342	0.9932	0.9986	74	3	2.5149	7.5447	−4.6388	0.9246	0.9704
	0.2	74	4	2.6845	10.7380	−3.8285	0.9945	0.9984	74	3	2.4851	7.4553	−4.8361	0.9262	0.9701
	0.3	74	5	2.6693	13.3467	−4.0233	0.9992	0.9969	74	3	2.4346	7.3038	−5.0338	0.9295	0.9689
	0.4	74	5	2.6839	13.4197	−4.2183	0.9994	0.9968	74	3	3.0484	9.1453	−5.2371	0.9920	0.9361
	0.5	74	5	2.6982	13.4909	−4.4134	0.9996	0.9967	74	3	3.0973	9.2918	−5.4455	0.9972	0.9323
	0.6	74	5	2.7121	13.5604	−4.6084	0.9998	0.9965	74	3	3.1230	9.3690	−5.6545	0.9989	0.9304
	0.7	74	5	2.7257	13.6285	−4.8036	0.9999	0.9963	74	3	3.1392	9.4175	−5.8638	0.9996	0.9292
	0.8	74	5	2.7390	13.6950	−4.9987	0.9999	0.9960	74	3	3.1503	9.4510	−6.0732	0.9999	0.9284
	0.9	74	5	2.7520	13.7600	−5.1939	1.0000	0.9958	74	3	3.1585	9.4755	−6.2827	1.0000	0.9279
	1.0	74	5	2.7648	13.8238	−5.3891	1.0000	0.9955	74	3	3.1648	9.4943	−6.4923	1.0000	0.9275

. The designs with $\kappa =0$ and $\kappa =1$ correspond to locally c- and D-optimal designs, respectively. From the table, we find that the optimal number of units remains constant, which is determined by the cost parameters, as we will discuss later. The results in Table 3 show us that $cD$-optimal designs depend on the planning values of $\mathbf{\theta }$, the weight $\kappa$, and the removal proportion p. With the increase of the removal proportion p, the number of inspections will decrease, but the length of the inspection interval will increase. The duration of the experiment will increase when the weight approaches 1. In addition, efficiencies of the $cD$-optimal designs with different $\kappa$ to the locally D- and c-optimal designs imply that the $cD$-optimal designs can perform very well in most cases. However, in some special cases, such as $p=0.3$ and $\mu =\widehat{\mu }-2s_{\widehat{\mu }},\sigma =\widehat{\sigma }-2s_{\widehat{\sigma }}$, the $cD$-optimal compound designs are very sensitive to the change of the weight parameter $\kappa$. Furthermore, with the increase of $\kappa$, the optimal number of inspections will increase or remain constant, and the duration for the design with $\kappa =1$ is always larger than that of the design with $\kappa =0$. To choose a proper value of the weight $\kappa$, we suggest using the efficiency lines plot. Figure 1 gives an illustrating example where efficiencies of the $cD$-optimal compound designs for different $\kappa$ to the locally D- and c-optimal designs are calculated and presented in a triangle or circle. Then, the horizontal coordinate of the intersection point of these two lines may be the best choice of $\kappa$.

In

Table 4. $BcD$-optimal PIC-I test plans for different weights $\kappa$ and different removal proportions p when the sets of the planning values are $\mathsf{\Theta }={\mathsf{\Theta }}_1$ and $\mathsf{\Theta }={\mathsf{\Theta }}_2$, respectively.

p	$\mathit{\kappa }$	${\mathbf{\Theta }}_1$							${\mathbf{\Theta }}_2$
		$\mathit{N}$	$\mathit{k}$	$\tau$	$\mathit{k}\tau$	${\mathsf{\Psi }}_{\mathit{BcD}}$	${\mathrm{Eff}}_{\mathit{BD}}\left({\xi }_{\mathit{BcD}}^{\ast }\right)$	${\mathrm{Eff}}_{\mathit{Bc}}\left({\xi }_{\mathit{BcD}}^{\ast }\right)$	$\mathit{N}$	$\mathit{k}$	$\tau$	$\mathit{k}\tau$	${\mathsf{\Psi }}_{\mathit{BcD}}$	${\mathrm{Eff}}_{\mathit{BD}}\left({\xi }_{\mathit{BcD}}^{\ast }\right)$	${\mathrm{Eff}}_{\mathit{Bc}}\left({\xi }_{\mathit{BcD}}^{\ast }\right)$
0.1	0.0	74	5	2.1731	10.8654	−3.5759	0.9360	1.0000	74	6	2.0791	12.4747	−3.6117	0.9673	1.0000
	0.1	74	7	2.1752	15.2266	−3.7842	0.9942	0.9977	74	7	2.1093	14.7654	−3.8200	0.9894	0.9989
	0.2	74	7	2.2187	15.5306	−3.9950	0.9959	0.9974	74	8	2.1187	16.9493	−4.0297	0.9966	0.9976
	0.3	74	7	2.2559	15.7910	−4.2059	0.9970	0.9970	74	8	2.1503	17.2027	−4.2399	0.9975	0.9973
	0.4	74	7	2.2885	16.0192	−4.4169	0.9978	0.9966	74	8	2.1789	17.4312	−4.4503	0.9981	0.9969
	0.5	74	8	2.2654	18.1234	−4.6283	0.9994	0.9955	74	8	2.2048	17.6381	−4.6607	0.9986	0.9965
	0.6	74	8	2.2880	18.3039	−4.8397	0.9996	0.9952	74	9	2.1921	19.7287	−4.8714	0.9996	0.9953
	0.7	74	8	2.3091	18.4726	−5.0511	0.9998	0.9949	74	9	2.2116	19.9047	−5.0821	0.9998	0.9950
	0.8	74	8	2.3289	18.6311	−5.2626	0.9999	0.9945	74	9	2.2299	20.0691	−5.2929	0.9999	0.9947
	0.9	74	8	2.3476	18.7806	−5.4742	1.0000	0.9942	74	9	2.2470	20.2229	−5.5038	1.0000	0.9944
	1.0	74	8	2.3653	18.9221	−5.6858	1.0000	0.9938	74	9	2.2630	20.3672	−5.7146	1.0000	0.9941
0.3	0.0	74	3	3.3461	10.0384	−3.4621	0.9323	1.0000	74	4	3.1924	12.7695	−3.4840	0.9817	1.0000
	0.1	74	4	3.3257	13.3029	−3.6553	0.9906	0.9994	74	5	3.2009	16.0044	−3.6769	0.9960	0.9987
	0.2	74	5	3.3230	16.6148	−3.8494	0.9986	0.9976	74	5	3.2300	16.1501	−3.8712	0.9969	0.9986
	0.3	74	5	3.3438	16.7190	−4.0442	0.9990	0.9974	74	5	3.2584	16.2922	−4.0656	0.9976	0.9983
	0.4	74	5	3.3641	16.8204	−4.2392	0.9993	0.9973	74	5	3.2860	16.4299	−4.2601	0.9982	0.9980
	0.5	74	5	3.3838	16.9191	−4.4341	0.9995	0.9971	74	5	3.3126	16.5628	−4.4547	0.9987	0.9976
	0.6	74	5	3.4031	17.0154	−4.6291	0.9997	0.9969	74	5	3.3381	16.6903	−4.6494	0.9990	0.9972
	0.7	74	5	3.4219	17.1094	−4.8242	0.9998	0.9966	74	5	3.3624	16.8121	−4.8441	0.9993	0.9967
	0.8	74	5	3.4403	17.2013	−5.0193	0.9999	0.9963	74	6	3.3533	20.1196	−5.0391	0.9999	0.9950
	0.9	74	5	3.4582	17.2911	−5.2144	1.0000	0.9960	74	6	3.3739	20.2431	−5.2341	1.0000	0.9946
	1.0	74	5	3.4758	17.3789	−5.4096	1.0000	0.9957	74	6	3.3936	20.3614	−5.4292	1.0000	0.9941

, we show $BcD$-optimal PIC-I test plans for different weights and removal proportions when the sets of the planning values are ${\mathsf{\Theta }}_1=[\widehat{\mu }-s_{\widehat{\mu }},\widehat{\mu }+s_{\widehat{\mu }}]\times [\widehat{\sigma }-s_{\widehat{\sigma }},\widehat{\sigma }+s_{\widehat{\sigma }}]$ and ${\mathsf{\Theta }}_2=[\widehat{\mu }-2s_{\widehat{\mu }},\widehat{\mu }+2s_{\widehat{\mu }}]\times [\widehat{\sigma }-2s_{\widehat{\sigma }},\widehat{\sigma }+2s_{\widehat{\sigma }}]$, respectively. Comparing the designs with $p=0.1$ and $p=0.3$, we find that when the removal proportion increases, the number of inspections will decrease and the lengths of the inspection time intervals will increase. Designs with $\mathsf{\Theta }={\mathsf{\Theta }}_1$ and $\mathsf{\Theta }={\mathsf{\Theta }}_2$ indicate that when the uncertainty of the planning values increase, the numbers of inspections tend to increase or remain constant, the lengths of the inspection intervals tend to decrease, and the durations of the optimal designs will become longer for most cases. Furthermore, we also compute the efficiencies of the $BcD$-optimal designs relative to the $BD$-optimal designs ($\kappa =1$) and to the $Bc$-optimal designs ($\kappa =0$) by the formula given in (13), in which $BD$- and $Bc$-optimality criteria are considered, respectively and then presented them in columns 8 (or 15) and 9 (or 16), respectively. By these efficiencies, the robust weights can be determined by the plot of efficiency lines, similar as we have shown in the $cD$-optimal compound designs.

5.3. Influence of the Cost Parameters

In order to investigate the influence of the cost parameters on the final robust optimal PIC-I test plans, we change the cost parameters over the sets $C_r\in \{4000,6000,8000\}$, $C_s\in \{70,80,90\}$, $C_i\in \{2,3,5\}$, $C_o\in \{1,2.5,4\}$ and compute the optimal test plans using Algorithm 1. The results for $BcD$-optimal designs with $p=0.3$ are presented in

Table 5. $BcD$-optimal PIC-I test plans for different cost parameters $C,C_s,C_i,C_o$ and weight parameters $\kappa$/$\omega$ when the removal proportion is $p=0.3$ and $\mathsf{\Theta }={\mathsf{\Theta }}_1$.

$\mathit{\kappa }$	C	${\mathit{C}}_{\mathit{s}}$	${\mathit{C}}_{\mathit{i}}$	${\mathit{C}}_{\mathit{o}}$	N	k	$\mathit{\tau }$	$\mathit{k}\mathit{\tau }$	${\mathbf{\Psi }}_{\mathit{BcD}}$	${\mathbf{Eff}}_{\mathit{D}}\left({\mathit{\xi }}_{\mathit{BcD}}^{\ast }\right)/$	${\mathbf{Eff}}_{\mathit{c}}\left({\mathit{\xi }}_{\mathit{BcD}}^{\ast }\right)$
0	4000	80	3	2.5	49	3	3.3415	10.0244	−3.0538	0.9924	0.9962
	6000				74	3	3.3461	10.0384	−3.4621	0.9917	0.9963
	8000				99	3	3.3485	10.0454	−3.7512	0.9913	0.9963
	6000	70	3	2.5	85	3	3.3461	10.0384	−3.5956	0.9917	0.9963
		80			74	3	3.3461	10.0384	−3.4621	0.9917	0.9963
		90			66	3	3.3461	10.0384	−3.3443	0.9917	0.9963
	6000	80	2	2.5	74	3	3.3461	10.0384	−3.4626	0.9915	0.9963
			3		74	3	3.3461	10.0384	−3.4621	0.9917	0.9963
			5		74	3	3.3461	10.0384	−3.4611	0.9920	0.9963
	6000	80	3	1.0	74	4	3.2950	13.1801	−3.4649	0.9956	0.9966
				2.5	74	3	3.3461	10.0384	−3.4621	0.9917	0.9963
				4.0	74	3	3.3405	10.0216	−3.4596	0.9923	0.9962
0.5	4000	80	3	2.5	49	5	3.3742	16.8709	−4.0238	0.9892	0.9896
	6000				74	5	3.3838	16.9191	−4.4341	0.9900	0.9914
	8000				99	5	3.3886	16.9432	−4.7242	0.9904	0.9923
	6000	70	3	2.5	84	5	3.3838	16.9191	−4.5677	0.9900	0.9914
		80			74	5	3.3838	16.9191	−4.4341	0.9900	0.9914
		90			66	5	3.3838	16.9191	−4.3164	0.9900	0.9914
	6000	80	2	2.5	74	5	3.3838	16.9192	−4.4350	0.9902	0.9917
			3		74	5	3.3838	16.9191	−4.4341	0.9900	0.9914
			5		74	5	3.3838	16.9190	−4.4325	0.9897	0.9907
	6000	80	3	1.0	74	5	3.3953	16.9767	−4.4384	0.9907	0.9930
				2.5	74	5	3.3838	16.9191	−4.4341	0.9900	0.9914
				4.0	73	5	3.3723	16.8615	−4.4299	0.9893	0.9898
1	4000	80	3	2.5	49	5	3.4664	17.3321	−4.9992	0.9847	0.9864
	6000				74	5	3.4758	17.3789	−5.4096	0.9856	0.9883
	8000				99	5	3.4805	17.4023	−5.6997	0.9860	0.9892
	6000	70	3	2.5	84	5	3.4758	17.3789	−5.5431	0.9856	0.9883
		80			74	5	3.4758	17.3789	−5.4096	0.9856	0.9883
		90			66	5	3.4758	17.3789	−5.2918	0.9856	0.9883
	6000	80	2	2.5	74	5	3.4758	17.3790	−5.4104	0.9858	0.9886
			3		74	5	3.4758	17.3789	−5.4096	0.9856	0.9883
			5		74	5	3.4758	17.3788	−5.4079	0.9853	0.9876
	6000	80	3	1.0	74	5	3.4870	17.4348	−5.4140	0.9863	0.9899
				2.5	74	5	3.4758	17.3789	−5.4096	0.9856	0.9883
				4.0	73	5	3.4646	17.3229	−5.4052	0.9849	0.9867

. The ranges of the planning values are $\mathsf{\Theta }={\mathsf{\Theta }}_1$. It is observed from Table 5 that the number of units will increase with an increase of the budget C, but decrease with the increase of the cost of one test unit $C_s$. The test plans are insensitive to the values of the cost parameter of one inspection $C_i$. However, the number of inspections or the duration of the experiment will decrease or remain constant with the increase of the cost operation $C_o$. Under different cost parameters, we also list the efficiencies of the $BcD$-optimal designs to the locally optimal design with the planning values $\mu =log5$ and $\sigma =0.5$. These efficiencies indicate that the performances of the robust $BcD$-optimal designs are actually quite insensitive to the cost parameters.

5.4. Optimal General PIC-I Test Plans

In this subsection, we consider the solution to the general problem given in (22). The model and values of the model parameters are the same as those described in the previous subsections. We use the PSO algorithm (Algorithm 2) to compute all general PIC-I test plans. Noting from the results in the previous tables that lengths of the inspection intervals are no more than 5, we then limit the lengths of the inspection intervals to the range of 0 to 10 in the PSO algorithm. The optimal PIC-I test plans are calculated for different design criteria and removal proportions. During the computation, we start Algorithm 2 by first setting $k=2$, calculate the k-point optimal inspection time points $t_i$ and then the optimal size of test units N, and increase k by 1 until the optimality criterion function does not decrease or arrives at the maximum inspection time $k_{max}$, see Algorithm 2 for details. We use the solutions obtained using the L-BFGS-B method as the initial values of the PSO algorithm. During the computation, we find that the design criterion is a convex function of k. See Figure 2 for an illustration. For the sake of saving space, we provide only the results with $p=0.3$.

Table 6. Locally D- and c-optimal general PIC-I test plans for different numbers of inspections with $\mu =log5$, $\sigma =0.5$ and $p=0.3$.

	k	N	Inspection Times	${\mathbf{\Psi }}_{\mathit{D}}/{\mathbf{\Psi }}_{\mathit{c}}$
	2	74	(2.612, 8.067)	−5.3216
	3	74	(2.547, 7.214, 9.929)	−5.4120
D-opt	4	74	(2.521, 6.961, 9.235, 11.333)	−5.4308
	5	74	(2.514, 6.896, 9.066, 10.767, 12.455)	−5.4344
	6	74	(2.512, 6.881, 9.029, 10.644, 12.002, 13.207)	−5.4346
	7	74	(2.512, 6.879, 9.025, 10.629, 11.948, 12.999, 13.434)	−5.4341
	2	74	(3.154, 5.554)	−3.4366
c-opt	3	74	(2.970, 4.886, 6.724)	−3.4581
	4	74	(2.934, 4.761, 6.216, 7.603)	−3.4610
	5	74	(2.928, 4.741, 6.140, 7.252, 8.083)	−3.4609

shows the locally D- and c-optimal PIC-I test plans for different number of inspections k and the values of their corresponding design criteria.

Table 7. $BcD$-optimal general PIC-I test plans for different numbers of inspections with $p=0.3$ and $\mathsf{\Theta }={\mathsf{\Theta }}_i,i=1,2$.

$\mathbf{\Theta }$	$\mathit{\kappa }$	k	N	Inspection Times	${\mathbf{\Psi }}_{\mathit{BcD}}$	${\mathbf{Eff}}_{\mathit{D}}\left({\mathit{\xi }}_{\mathit{BcD}}^{\ast }\right)$	${\mathbf{Eff}}_{\mathit{c}}\left({\mathbf{\xi }}_{\mathit{BcD}}^{\ast }\right)$
${\mathsf{\Theta }}_1$	0	2	74	(3.911, 7.256)	−3.4467	0.8388	0.8882
		3	74	(3.700, 6.219, 8.875)	−3.4764	0.8693	0.9238
		4	74	(3.658, 6.020, 8.051, 10.269)	−3.4811	0.8791	0.9314
		5	74	(3.650, 5.983, 7.900, 9.628, 11.190)	−3.4813	0.8818	0.9329
		6	74	(3.650, 5.982, 7.892, 9.595, 11.047, 11.312)	−3.4808	0.8816	0.9329
	0.5	2	74	(3.811, 9.328)	−4.3429	0.8983	0.8882
		3	74	(3.540, 7.755, 11.515)	−4.4152	0.9313	0.9088
		4	74	(3.453, 7.255, 10.413, 13.462)	−4.4335	0.9383	0.9218
		5	74	(3.430, 7.122, 10.123, 12.563, 15.151)	−4.4376	0.9407	0.9256
		6	74	(3.424, 7.090, 10.053, 12.353, 14.409, 16.276)	−4.4378	0.9414	0.9263
		7	74	(3.423, 7.087, 10.046, 12.334, 14.342, 16.028, 16.511)	−4.4373	0.9415	0.9263
	1	2	74	(3.304, 9.810)	−5.3200	0.8949	0.8700
		3	74	(3.227, 8.798, 12.277)	−5.4196	0.9254	0.8763
		4	74	(3.193, 8.480, 11.343, 14.299)	−5.4422	0.9367	0.8813
		5	74	(3.183, 8.389, 11.098, 13.437, 16.036)	−5.4472	0.9405	0.8829
		6	74	(3.180, 8.366, 11.037, 13.234, 15.279, 17.240)	−5.4476	0.9414	0.8831
		7	74	(3.180, 8.363, 11.030, 13.209, 15.188, 16.899, 17.583)	−5.4472	0.9415	0.8831
${\mathsf{\Theta }}_2$	0	2	74	(3.835, 7.571)	−3.4431	0.8794	0.8994
		3	74	(3.590, 6.233, 9.451)	−3.4909	0.8914	0.9388
		4	74	(3.538, 5.963, 8.319, 11.210)	−3.4998	0.8963	0.9472
		5	74	(3.527, 5.902, 8.074, 10.268, 12.605)	−3.5010	0.8986	0.9490
		6	74	(3.525 5.892 8.035 10.126 12.055 13.172)	−3.5006	0.8993	0.9492
	0.5	2	74	(3.713, 8.514)	−4.3305	0.9316	0.9037
		3	74	(3.469, 7.312, 11.248)	−4.4273	0.9320	0.9233
		4	74	(3.389, 6.910, 10.013, 13.540)	−4.4508	0.9390	0.9341
		5	74	(3.364, 6.784, 9.666, 12.385, 15.614)	−4.4566	0.9424	0.9377
		6	74	(3.357, 6.750, 9.572, 12.096, 14.599, 17.120)	−4.4573	0.9435	0.9384
		7	74	(3.356,6.745,9.559,12.055,14.459,16.629,17.601)	−4.4569	0.9436	0.9384
	1	2	74	(3.235, 8.905)	−5.2680	0.9631	0.8940
		3	74	(3.157, 7.951, 11.858)	−5.4075	0.9569	0.9023
		4	74	(3.120, 7.644, 10.724, 14.294)	−5.4409	0.9634	0.9082
		5	74	(3.106, 7.544, 10.403, 13.144, 16.511)	−5.4492	0.9667	0.9103
		6	74	(3.101, 7.515, 10.315, 12.851, 15.443, 18.249)	−5.4506	0.9677	0.9106
		7	74	(3.101, 7.510, 10.298, 12.796, 15.253, 17.559, 19.001)	−5.4503	0.9679	0.9105

shows the $BcD$-optimal PIC-I test plans with different $\kappa$ when the regions of parameters are $\mathsf{\Theta }={\mathsf{\Theta }}_1$ and $\mathsf{\Theta }={\mathsf{\Theta }}_2$, respectively. In Table 6 and Table 7, the optimal test plans are shown in bold. Looking at these tables we can see that the number of time points of every optimal PIC-I test plan does not exceed 10, which is very similar to its corresponding equispaced design. In Table 7, the efficiencies of k points $BcD$-optimal designs relative to their corresponding locally ones given in Table 6 are also given in the last two columns. The results indicate that the $BcD$-optimal test plans can perform very well in most situations, especially when the planning values are in the prior regions. Finally, to investigate the influence of the cost parameters on the optimal test plans, the $BcD$-optimal designs are also calculated under different combinations of the cost parameters, and the results for $p=0.3$ are presented in

Table 8. $BcD$-optimal general PIC-I test plans for different cost parameters $C,C_s,C_i,C_o$ and weight parameters $\kappa$ when the removal proportion $p=0.3$ and $\mathsf{\Theta }={\mathsf{\Theta }}_1$.

$\mathit{\kappa }$	C	${\mathit{C}}_{\mathit{s}}$	${\mathit{C}}_{\mathit{i}}$	${\mathit{C}}_{\mathit{o}}$	N	Inspection Times	${\mathbf{\Psi }}_{\mathit{BcD}}$
	4000	80	3	2.5	49	(3.659, 6.020, 8.035, 10.159)	−3.0725
	6000				74	(3.650, 5.983, 7.900, 9.628, 11.190)	−3.4813
	8000				99	(3.650, 5.982, 7.902, 9.658, 11.371)	−3.7708
	6000	70	3	2.5	85	(3.650, 5.983, 7.900, 9.628, 11.190)	−3.6148
		80			74	(3.650, 5.983, 7.900, 9.628, 11.190)	−3.4813
0		90			66	(3.650, 5.983, 7.900, 9.628, 11.190)	−3.3635
	6000	80	2	2.5	74	(3.650, 5.983, 7.900, 9.628, 11.190)	−3.4821
			3		74	(3.650, 5.983, 7.900, 9.628, 11.190)	−3.4813
			5		74	(3.658, 6.020, 8.051, 10.269)	−3.4798
	6000	80	3	1.0	74	(3.649, 5.981, 7.908, 9.712, 11.677)	−3.4842
				2.5	74	(3.650, 5.983, 7.900, 9.628, 11.190)	−3.4813
				4.0	74	(3.659, 6.020, 8.032, 10.138)	−3.4785
	4000	80	3	2.5	49	(3.428, 7.115, 10.106, 12.509, 14.948)	−4.0277
	6000				74	(3.424, 7.090, 10.053, 12.353, 14.409, 16.276)	−4.4378
	8000				99	(3.425, 7.092, 10.058, 12.372, 14.478, 16.538)	−4.7280
	6000	70	3	2.5	84	(3.424, 7.090, 10.053, 12.353, 14.409, 16.276)	−4.5714
		80			74	(3.424, 7.090, 10.053, 12.353, 14.409, 16.276)	−4.4378
0.5		90			66	(3.424, 7.090, 10.053, 12.353, 14.409, 16.276)	−4.3200
	6000	80	2	2.5	74	(3.424, 7.090, 10.053, 12.353, 14.409, 16.277)	−4.4388
			3		74	(3.424, 7.090, 10.053, 12.353, 14.409, 16.276)	−4.4378
			5		74	(3.430, 7.122, 10.123, 12.562, 15.151)	−4.4359
	6000	80	3	1.0	74	(3.426, 7.096, 10.068, 12.403, 14.592, 16.976)	−4.4420
				2.5	74	(3.424, 7.090, 10.053, 12.353, 14.409, 16.276)	−4.4378
				4.0	73	(3.422, 7.085, 10.042, 12.316, 14.271, 15.755)	−4.4338
	4000	80	3	2.5	49	(3.182, 8.384, 11.085, 13.392, 15.858)	−5.0371
	6000				74	(3.180, 8.366, 11.037, 13.234, 15.279, 17.240)	−5.4476
	8000				99	(3.181, 8.367, 11.042, 13.250, 15.341, 17.480)	−5.7379
	6000	70	3	2.5	84	(3.180, 8.366, 11.037, 13.234, 15.279, 17.240)	−5.5812
		80			74	(3.180, 8.366, 11.037, 13.234, 15.279, 17.240)	−5.4476
		90			65	(3.180, 8.366, 11.037, 13.234, 15.279, 17.240)	−5.3299
1	6000	80	2	2.5	74	(3.180, 8.366, 11.037, 13.234, 15.279, 17.241)	−5.4487
			3		74	(3.180, 8.366, 11.037, 13.234, 15.279, 17.240)	−5.4476
			5		74	(3.180, 8.366, 11.037, 13.234, 15.279, 17.238)	−5.4456
	6000	80	3	1.0	74	(3.181, 8.370, 11.049, 13.277, 15.442, 17.875)	−5.4521
				2.5	74	(3.180, 8.366, 11.037, 13.234, 15.279, 17.240)	−5.4476
				4.0	73	(3.179, 8.363, 11.028, 13.201, 15.153, 16.760)	−5.4433

. Similar conclusions to their corresponding ES designs can be drawn from the table.

Algorithm 2 PSO Algorithm.

1: Set the values of the cost parameters $C_s,C_i,C_o,C$; design space related parameters $k^{\diamond },p_j$; the set of values of the parameters $\mathbf{\beta }$, $\mathsf{\Theta }$; the hyperparamters in the priors for Bayesian optimality criterion; the weight parameter $\kappa$ or $\omega$. 2: Set $k=2$. 3: Find the k-point optimal design ${\xi }_k^{\ast }$ by PSO. 3.1: Initialize particles. Initialize n random particles $X_i^k$ and velocities $V_i^k,i=1,\dots ,n$. Calculate the fitness values $f\left(X_i^k\right)$. Determine local and global best positions $P_i^k$ and $P_g^k$. 3.2: Repeat until the stopping criterion is met. 3.2.1: Calculate particle velocity according the following equation $$V_i^k=\alpha V_i^k+{\alpha }_1{u}_1(P_i^k-X_i^k)+{\alpha }_2{u}_2(P_g^k-X_i^k).$$ 3.2.2: Update the variable vector $X_i^k$ according the following equation $$X_i^k=X_i^k+V_i^k.$$ 3.2.3: Calculate fitness values $f\left(X_i^k\right)$. 3.2.4: Update best positions $P_i^k$ and $P_g^k$. 3.3: Output ${\xi }_k^{\ast }=argminf\left(X^k\right)$ with $gbes{t}_k=f\left({\xi }_k^{\ast }\right)$. 4: Set $k=k+1$. If $k\le k^{\diamond }$ or if the other stop rule is not satisfied, go to Step 3, else go to Step 5. 5: Output ${\xi }^{\ast }=\{N^{\ast },arg{min}_{k}f\left({\xi }_k^{\ast }\right)$} with $gbest=f\left({\xi }^{\ast }\right)$.

6. A Real Life Example

In this section we consider a real data example, which contains 112 patients with plasma cell myeloma treated at the National Cancer Institute (Carbone et al. [37]) to illustrate our method. For easy reference, we reproduce this data set in

Table 9. Plasma cell myeloma survival times.

Interval in Months	Number at Risk	Number of Withdrawals
$[0,5.5)$	112	1
$[5.5,10.5)$	93	1
$[10.5,15.5)$	76	3
$[15.5,20.5)$	55	0
$[20.5,25.5)$	45	0
$[25.5,30.5)$	34	1
$[30.5,40.5)$	25	2
$[40.5,50.5)$	10	3
$[50.5,60.5)$	3	2
$[60.5,\infty )$	0	0

. Note that this data set has been reanalysed by many authors (see Ng and Wang [6] and Lin and lio [8]) via different distributions and estimating methods (Balakrishnan and Cramer [2], Ch. 18). To be consistent with the topic we investigated in this paper, we fit the data by the Weibull distribution (1), similar as Ng and Wang [6], and the parameters are estimated by the maximum likelihood method. The resulted estimates are given as follows: $\widehat{\mu }=3.1391,\widehat{\sigma }=0.8132$ and $s_{\widehat{\mu }}=0.0841,s_{\widehat{\sigma }}=0.0724$. Similar in the simulation example, we assume that the range of the values of $\mu$ and $\sigma$ are $[a,b]=[\widehat{\mu }-s_{\widehat{\mu }},\widehat{\mu }+s_{\widehat{\mu }}]$ and $[c,d]=[\widehat{\sigma }-s_{\widehat{\sigma }},\widehat{\sigma }+s_{\widehat{\sigma }}]$, respectively. Furthermore, we assume that the hyperparameters in the prior distributions $N({\mu }_0,{\sigma }_0^2)$ and ${\mathsf{\Gamma }}^{-1}({\nu }_0,{\gamma }_0)$ are ${\mu }_0=\widehat{\mu },{\sigma }_0=\widehat{\sigma },{\nu }_0={\widehat{\sigma }}^2/s_{\widehat{\sigma }}^2+2=128.1588,{\gamma }_0=\widehat{\sigma }({\widehat{\sigma }}^2/s_{\widehat{\sigma }}^2+1)=103.4055$. Since we have no any cost information for the original experiment, we then assume that $C_s=\$60$/unit, $C_i=\$30$, $C_o=\$25$/m. Therefore, the total cost is $C=121C_s+9C_i+60.5C_0=\$9042.5$. In $BcD$- or $McD$-optimal design criteria, we set $q=0.5$, which means that the median lifetime of survivors is of interest. In order to compare with the original inspection scheme ${\xi }_{\mathrm{org}}$, we first calculate the optimal ES and the general PIC-I test plans, respectively, assuming that the number of inspection times and the removal scheme are the same as in ${\xi }_{\mathrm{org}}$. The results for different design criteria are shown in

Table 10. Optimal ES PIC-I test plans in the real data example when the removal scheme is the same as that in the original design.

Design Criterion	$\mathit{\kappa }/\mathit{\omega }$	N	$\mathit{\tau }$	$\mathit{k}\mathit{\tau }$	${\mathbf{\Psi }}_{\mathbf{L}}$	${\mathbf{Eff}}_{\mathbf{L}}\left({\mathit{\xi }}_{\mathbf{org}}\right)$
D-opt	-	132	3.5972	32.3744	−9.6813	0.8243
c-opt	-	134	3.1923	28.7311	−9.6071	0.8233
BcD-opt	0.0	128	4.8029	43.2262	−5.0343	1.0001
	0.5	123	5.9594	53.6349	−5.1176	1.0086
	1.0	120	6.7642	60.8781	−5.2223	0.9957
McD-opt	0.0	128	4.7736	42.9627	−4.8030	1.0061
	0.5	122	6.2919	56.6271	−4.8794	1.0132
	1.0	118	7.3139	65.8255	−4.9846	0.9914

and

Table 11. Optimal general PIC-I test plans in the real data example when the removal scheme is the same as that in the original design.

Design Criterion	$\mathit{\kappa }/\mathit{\omega }$	N	$\mathit{\xi }$	${\mathbf{\Psi }}_{\mathbf{L}}$	${\mathbf{Eff}}_{\mathbf{L}}\left({\mathit{\xi }}_{\mathbf{org}}\right)$
D-opt	-	134	(17.524, 19.856, 21.544, 22.897, 24.076, 25.076, 26.076, 27.076, 28.076)	−9.9818	0.6103
c-opt	-	134	(16.378, 18.478, 19.831, 20.831, 21.831, 22.831, 25.458, 26.458, 27.458)	−9.8213	0.6646
BcD-opt	0.0	124	(1.089, 3.549, 7.295, 11.990, 19.281, 39.281, 52.830, 52.930, 53.030)	−5.1580	0.8838
	0.5	122	(1.239, 4.075, 8.557, 14.657, 25.349, 42.696, 57.013, 57.113, 57.213)	−5.2252	0.9058
	1.0	120	(1.568, 5.263, 11.588, 21.011, 34.291, 47.815, 60.597, 60.697, 60.797)	−5.3025	0.9190
McD-opt	0.0	122	(0.856, 3.094, 6.782, 11.590, 19.249, 39.249, 55.940, 56.040, 56.140)	−4.9486	0.8698
	0.5	120	(0.981, 3.582, 8.018, 14.255, 25.322, 45.322, 62.099, 62.199, 62.299)	−5.0123	0.8871
	1.0	118	(1.250, 4.668, 10.974, 20.715, 36.007, 52.237, 66.919, 67.019, 67.119)	−5.0858	0.8960

, respectively. Efficiencies of ${\xi }_{\mathrm{org}}$ under different design criteria are also given in the last column of these tables. These results indicate that the original inspection scheme ${\xi }_{\mathrm{org}}$ has good robustness under different design criteria, except general D- and c-optimal design criteria. In addition, the optimal general PIC-I test plan has better performance than the corresponding ES test plan. Furthermore, considering different removal proportions, $p=0$ and $p=0.3$, we redesign the Plasma cell myeloma inspection scheme. Optimal robust ES and general inspection schemes for different design criteria are shown in

Table 12. Optimal ES PIC-I test plans in the real data example.

	p	$\mathit{\kappa }$	N	k	$\mathit{\tau }$	$\mathit{k}\mathit{\tau }$	${\mathbf{\Psi }}_{\mathbf{L}}$	${\mathbf{Eff}}_{\mathbf{L}}\left({\mathit{\xi }}_{\mathbf{org}}\right)$
D-opt	0	-	127	24	1.1200	26.8790	−9.9947	0.6025
	0.3	-	138	4	6.2843	25.1371	−8.6852	2.2319
c-opt	0	-	129	21	1.2749	26.7728	−9.8208	0.6649
	0.3	-	137	3	9.3286	27.9858	−8.8469	1.7608
BcD-opt	0	0.0	120	19	2.5650	48.7345	−5.1274	0.9112
		0.5	119	17	3.2444	55.1550	−5.2181	0.9122
		1.0	117	16	3.7316	59.7058	−5.3183	0.9046
	0.3	0.0	134	7	4.4587	31.2107	−4.7950	1.2705
		0.5	127	6	8.2476	49.4854	−4.7230	1.4966
		1.0	126	4	13.1284	52.5135	−4.7846	1.5425
McD-opt	0	0.0	120	19	2.6290	49.9512	−4.9044	0.9091
		0.5	117	18	3.2852	59.1339	−4.9892	0.9079
		1.0	115	16	4.0608	64.9733	−5.0872	0.8947
	0.3	0.0	147	3	1.3227	3.9682	−4.5884	1.2470
		0.5	127	6	8.0017	48.0105	−4.4663	1.5315
		1.0	125	4	14.1928	56.7714	−4.5415	1.5441

and

Table 13. Optimal general PIC-I test plans in the real data example.

	p	$\mathit{\kappa }$	N	$\mathit{\xi }$	${\mathbf{\Psi }}_{\mathbf{L}}$	${\mathbf{Eff}}_{\mathbf{L}}\left({\mathit{\xi }}_{\mathbf{org}}\right)$
D-opt	0	-	133	(15.169, 17.295, 18.779, 19.956, 20.972, 21.972, 22.972, 23.972, 24.972, 25.972, 26.972)	−10.0496	0.5703
	0.3	-	137	(20.563, 24.384, 25.548, 26.548)	−9.6097	0.8854
c-opt	0	-	133	(14.391, 16.419, 17.831, 18.935, 19.935, 20.935, 21.935, 22.935, 25.091, 26.091, 27.091)	−9.8711	0.6323
	0.3	-	137	(19.446, 21.483, 25.585, 26.585)	−9.6065	0.8238
BcD-opt	0	0.0	123	(0.303, 1.086, 2.409, 4.305, 6.836, 10.150, 14.604, 21.274, 41.274, 53.615)	−5.2035	0.8444
		0.5	121	(0.400, 1.427, 3.154, 5.647, 9.058, 13.761, 20.764, 32.648, 45.421, 58.097)	−5.2822	0.8556
		1.0	120	(0.514, 1.826, 4.045, 7.325, 12.036, 19.010, 28.766, 39.155, 49.847, 61.535)	−5.3693	0.8596
	0.3	0.0	128	(4.020, 10.849, 30.849, 48.098)	−4.9285	1.1117
		0.5	126	(4.899, 16.739, 36.739, 52.615)	−4.8890	1.2677
		1.0	125	(6.771, 26.771, 46.771, 56.246)	−4.9230	1.3431
McD-opt	0	0.0	122	(0.214, 0.856, 2.035, 3.832, 6.349, 9.776, 14.528, 21.800, 41.800, 56.687)	−4.9936	0.8315
		0.5	119	(0.275, 1.095, 2.598, 4.902, 8.192, 12.857, 19.849, 31.909, 47.622, 63.043)	−5.0677	0.8393
		1.0	117	(0.361, 1.432, 3.407, 6.500, 11.122, 18.163, 28.776, 41.060, 53.851, 67.732)	−5.1514	0.8391
	0.3	0.0	127	(3.524, 10.215, 30.215, 50.215)	−4.7243	1.0189
		0.5	124	(5.391, 25.391, 45.391, 57.616)	−4.6739	1.2444
		1.0	123	(6.215, 26.215, 46.215, 60.490)	−4.6969	1.3219

, respectively. The efficiencies of the original scheme ${\xi }_{\mathrm{org}}$ are also listed in the last columns of these tables. These numerical results suggest that it is important to consider the removal proportions (or dropouts) at the design stage. The structures of robust PIC-I test schemes obtained under different design criteria may differ substantially.

7. Conclusions and Discussions

This paper considers robust optimal life testing plans under the PIC-I scheme. Based on the D-optimality and c-optimality, both the Bayesian and the minimax compound optimality criteria are proposed to obtain robust designs against the uncertainty of the model parameters and the design objectives. To make the designs more practical, the cost of the experiment at the design stage is taken into account. Two algorithms, including the MNO algorithm and the PSO algorithm, are proposed to solve the optimization problems. A numerical example and a real data example are provided to validate the feasibility and effectiveness of our methods. It is easy to see that the inspection times, the length of the inspection interval, and the experimental duration in the c- ($Bc$- or $Mc$-) optimal ES PIC-I test plan are less than the corresponding D- ($BD$- or $MD$-) optimal PIC-I test plan. Those in the $cD$- ($BcD$- or $McD$-) optimal PIC-I test plan are somewhere in between. Increasing the removal proportion will reduce the number of the inspections and lengthen the inspection interval. With the increase of the uncertainty of the model parameters, the number of the inspections and the duration of the optimal PIC-I test plan tend to increase. The $MCD$-optimal PIC-I test plan (not shown in the simulation example) tends to have more inspection times and longer duration of the experimental than its corresponding $BcD$-optimal test plan. Similar conclusions can be drawn for the general optimal test plan. In addition, we should point out that the optimal PIC-I test plans depend on the budget of the experiment. We need to carefully set the larger coefficients in the cost function, because they will have a big impact on the experiment. Overall, our numerical results indicate that the proposed approaches can obtain better designs with improved robustness for conducting PIC-I life tests to estimate the model parameters and the qth quantile of the lifetime distribution efficiently. They also provide a coherent way to efficiently use one’s resources. Our approaches are intuitive and can be useful to engineers.

In this paper, to approximate the Bayesian and minimax design criteria, we choose finite percentiles of the prior parameter distributions. To improve the accuracy of these approximations, some Monte Carlo sampling methods and advanced optimization methods (Ruidas et al. [33]) need to be applied. This may be a topic worth studying in the future. This paper mainly focuses on obtaining robust designs under PIC-I test schemes when the model parameters are unknown and there are more than one design objective. However, there are still other uncertainties that need to be considered in the design stage. For example, many studies on designing PIC-I test plans assume that the lifetime data follow some symmetrical or asymmetrical distribution, such as Weibull (Wu et al. [14]), lognormal (Roy and Pradhan [17,19]), truncated normal (Lodhi and Tripathi [10]), and Exponentiated Frech’et (Wu and Chang [12]). In fact, many symmetric and asymmetric distributions have recently been proposed to fit lifetime data and evaluate the reliability of system (see [38,39,40,41]). When life time data is available, many methods of hypothesis testing have been provided to check the fitness of a given distribution (Jäntschi [42]). However, we need to point out that the life data is not available at the design stage, so the latent lifetime distribution may be uncertain. Therefore, it is an interesting topic to find robust designs against possible departures from underlying model assumptions (Zhao et al. [43]). In addition, our methods proposed in this paper can be extended to the situations of accelerated life test (Wu and Huang [44]) or degradation life test [45].