**Integration of Dimension Reduction and Uncertainty Quantification in Designing Stretchable Strain Gauge Sensor**

#### **Sungkun Hwang 1, Recep M. Gorguluarslan 2, Hae-Jin Choi 3,\* and Seung-Kyum Choi 1,\***


Received: 29 November 2019; Accepted: 13 January 2020; Published: 16 January 2020

**Abstract:** Interests in strain gauge sensors employing stretchable patch antenna have escalated in the area of structural health monitoring, because the malleable sensor is sensitive to capturing strain variation in any shape of structure. However, owing to the narrow frequency bandwidth of the patch antenna, the operation quality of the strain sensor is not often assured under structural deformation, which creates unpredictable frequency shifts. Geometric properties of the stretchable antenna also severely regulate the performance of the sensor. Especially rugged substrate created by printing procedure and manual fabrication derives multivariate design variables. Such design variables intensify the computational burden and uncertainties that impede reliable analysis of the strain sensor. In this research, therefore, a framework is proposed not only to comprehensively capture the sensor's geometric design variables, but also to effectively reduce the multivariate dimensions. The geometric uncertainties are characterized based on the measurements from real specimens and a Gaussian copula is used to represent them with the correlations. A dimension reduction process with a clear decision criterion by entropy-based correlation coefficient dwindles uncertainties that inhibit precise system reliability assessment. After handling the uncertainties, an artificial neural network-based surrogate model predicts the system responses, and a probabilistic neural network derives a precise estimation of the variability of complicated system behavior. To elicit better performance of the stretchable antenna-based strain sensor, a shape optimization process is then executed by developing an optimal design of the strain sensor, which can resolve the issue of the frequency shift in the narrow bandwidth. Compared with the conventional rigid antenna-based strain sensors, the proposed design brings flexible shape adjustment that enables the resonance frequency to be maintained in reliable frequency bandwidth and antenna performance to be maximized under deformation. Hence, the efficacy of the proposed design framework that employs uncertainty characterization, dimension reduction, and machine learning-based behavior prediction is epitomized by the stretchable antenna-based strain sensor.

**Keywords:** stretchable antenna-based strain sensor; structural optimization; structural health monitoring; dimension reduction; entropy-based correlation coefficient; multidisciplinary design and analysis; uncertainty-integrated and machine learning-based surrogate modeling

#### **1. Introduction**

Structural health monitoring (SHM) is implemented to evaluate the physical conditions of structures with consistent surveillance. In accordance with the information garnered by an SHM system, engineers can identify the critical roots in structural damage or deterioration and provide applicable approaches to avoid structural failures. Compared with traditional fixed inspections with interval schedules that demand excessive maintenance [1], the SHM systems can operate for condition-based maintenance and reduce preventive maintenance cost, life cycle costs, and potential catastrophic failure [2]. Therefore, the systems can heighten both the capability and reliability of the monitoring system. Application of the SHM systems has been broadened from heavy mechanical equipment or civil structures to aerospace or bio-mechanical examples [3].

In the SHM systems, a strain is regarded as a vital factor, which should be rigorously examined. By measuring strain of structures, a strain gauge is able to predict certain failure modes such as crack propagation, deformation [4], vibration, or mechanical loading [5]. Between copious types of strain gauges, metallic foil pattern is mainly employed because it identifies electrical resistance connected to strain changes under structural deformation. In spite of the benefits of a foil gauge (e.g., simple circuit structure, lower fabrication cost, and applicability in various examples), its usability has been restrained owing to the reliance on long cables corresponding to the power supply and data transmission. An acceptable solution to the issue is to use sensors of passive wireless strain gauge, which reduces dependency on batteries and has lower installation costs. In recent years, great attention has been paid to passive wireless antennas on account of their compelling suitability to observe strain of structures [6].

An antenna installed on the surface of a deformable structure transfers details about the strain variations. It is approachable to evaluate the stability and physical circumstance of the structure with the transferred data of the strain variations. In the strain gauge markets, interests in micro-strip patch (MSP) antenna have increased because of lower fabrication cost, lighter load, lower profile planar configuration, and capability of multiband operation [7]. Thus, the majority have still exploited a rigid patch antenna to capture strain variation. However, the conventional rigid patch antenna often fails to deliver reliable information owing to two critical drawbacks: mediocre consideration for non-uniform geometric uncertainties [8] and limited bandwidth of resonance frequency (Rf) [9]. To assure the performance of the MSP antenna-based strain sensor, these two drawbacks should be meticulously examined.

The first drawback is about the geometric uncertainties on the surface of the substrate of an antenna shown in Figure 1a. As depicted in Figure 1b, a typical scanning electron microscope (SEM) shows that the substrate surface has certain scratched patterns, which makes the substrate rugged, so there exist obtrusive thickness variations on the substrate. Variation of the substrate roughness is unavoidable because the unstable fabrication process requires manual operation by engineers [10]. The geometric uncertainties due to this surface roughness are more critical for the stretchable antenna-based strain sensors, which broaden the applicability of such sensors for the wearable devices or human body. Compared with conventional rigid antennas, the uncertainties have a much stronger influence on the performance of the antenna along with the non-uniform changes under deformation [11–13]. Therefore, the substrate fabrication has been steadily operated with polydimethylsiloxane (PDMS) or other stretchable materials to be firmly attached to the surface of rough structures [14].

**Figure 1.** Stretchable antenna. (**a**) The antenna (dark color) and the polydimethylsiloxane (PDMS) substrate (light color); (**b**) scratched surface of the substrate observed by FEI Quanta 250 scanning electron microscope (SEM) device.

In the computer-based design and analysis of stretchable antennas for a strain gauge, usually, an assumption of uniform thickness is made to reduce design complexity and computational burden. However, if the randomness corresponding to the thickness is overlooked, reliable response assessment of the MSP antenna-based strain sensor can be impeded. As already proven by various simulations and measurements [15–20], although the analysis of the non-uniform substrate-based antenna is always more complicated than that of constant substrate-based antenna, it involves better reliability predictions as the approach corresponding to the uneven substrate represents the real environment. Thus, it is absolutely imperative to examine the thickness variation to improve the reliability of the antenna-based sensor.

Unfortunately, however, considering the roughness would bring extremely increased computational costs as well as total number of design variables (DVs). It is indispensable to capture and propagate all possible critical variables and the corresponding uncertainties to ensure that the system is reliable, but such an approach would incur high dimensions of DVs, increasing computational complexity [21]. Furthermore, once the DVs are contaminated by the immensely correlated random behavior, the extravagant variables obstruct the reliability of the performance estimation [22,23]. To eliminate these issues, a rigorous design of stretchable electronics based on efficient capturing and modeling of correlated and high dimensional random variables is required [24]. Therefore, this study proposed a dimension reduction (DR) framework to effectively manage uncertainties especially occurring for the substrate thickness. In the framework, a Gaussian copula is utilized [25] to precisely model the intricacy of DVs. The intricacy can be described by a joint distribution that consists of manifold marginal distributions. After modeling the geometric uncertainties, a clear guideline engaging the entropy-based correlation coefficient was exploited to suggest a better decision between feature extraction and selection. Through the guideline, engineers can embrace a competent choice of an efficient DR method based on multivariate data properties. In accordance with the entropy-based correlation coefficient, two DR methods, particularly feature extraction and feature selection, diminish the overabundance of the DVs. Specifically, the performance of two feature extraction methods (principal component analysis (PCA) [26] and auto-encoder (AE) [27]) are evaluated to identify the one that gives the best predictions, while feature selection employs independent features test (IFT) [28]. Details will be addressed in the following sections.

The second drawback is about the narrow bandwidth, which often fails to warrant the operation of the MSP antenna for measuring strain. In general, information transferred by a radio frequency system is only trustworthy within the desired frequency bandwidth. Owing to the unforeseen structural deformation that produces extreme frequency shifts, the stretchable strain sensor cannot perform stable detection when the frequency deviates from an allowable bandwidth. To handle this issue, manifold solutions have been proposed. Specifically, the ultra-wideband (UWB) wireless system is used to facilitate multiple broadband and a high data rate, which enlarges the range for frequency activity [29]. However, the solution is not convincing owing to the complicated design process, power limitation, and standardization. Other approaches engaging optimization techniques also exist, but most have focused on sizing adjustment of the antenna [30]. Even though the sizing optimization can draw simple and efficient antenna update [31], the technique forces the antenna to maintain its original shape, which interrupts a fundamental solution regarding the narrow bandwidth. Furthermore, topology optimization with solid isotropic material with penalization or the level set method has been exploited to address the issue [32], but the efficiency of the suggested design is only identified when it is fabricated by special additive manufacturing process, which intensely increases the fabrication burden. Assorted studies have also suggested unique shape modification such as conical, pentagon, or even fractal geometry [33,34], but such modifications refuse to preserve an initial outline of the MSP antenna that has advantages of cheaper and easier fabrication process and diverse applicability, so their efficiency is only meaningful in a specialized example.

To address these issues of the existing design approaches for the stretchable antenna, a design optimization process that includes a structural shape optimization is proposed in this study. In comparison with the existing approaches, the proposed design process not only inherits the advantages of the MSP antenna, but also derives an innovative antenna shape to enhance functionality in the multi-physics domain. Owing to the mechanical flexibility and electrical radio frequency behavior, the stretchable MSP antenna for strain gauge requires a meticulous multidisciplinary design optimization. Thus, the proposed design maximizes the frequency stability and antenna performance, referring to a return loss under unexpected structural deformation. The optimum design is garnered in accordance with the basic shape of the MSP antenna, so that it can truncate an intense fabrication burden. The design also operates at dual bandwidth to keep the benefit of UWB, but it induces less design intricacy. For the system response prediction, the proposed framework employs the artificial neuron network (ANN)-based surrogate model under geometric uncertainties stemming from the stretchability and fluctuating antenna substrate. Finally, the reliability of the proposed design is evaluated by the probabilistic neuron network (PNN) and Monte Carlo simulation (MCS).

This paper is structured as follows. Section 2 introduces backgrounds of the Rf fluctuation under structural deformation and the DR methods. Section 3 introduces the proposed design framework employing the DR methods and shape optimization to determine the optimal new design of the stretchable antenna-based strain sensor. In Section 4, an example of the MSP antenna employed as a stretchable strain gauge is provided to represent the efficacy and merits of the proposed framework.

#### **2. Micro-Strip Patch Antenna and DR Process**

#### *2.1. Micro-Strip Patch Antenna*

An MSP antenna is composed of a patch, a feed line, a dielectric substrate, and a ground, as depicted in Figure 2a. Owing to the effect of the inductance and the capacitance, intrinsic resonance frequency (Rf) is produced where inductive and capacitive reactance offset each other, which results in a transfer of signal and power. The Rf fluctuates as the antenna geometry is changed, which directly indicates how much the antenna is distorted. The antenna employed in a strain gauge normally adheres to the fragile surface of a structure, and the antenna's shape will be changed in which the structure is deformed. Rf is then promptly changed in accordance with the structural deformation of the antenna. Therefore, it is ineludible to monitor the Rf variation of the antenna to observe the status of the structured deformities [4]. As shown in Figure 2b, the Rf can be appraised based on the geometric information as follows:

$$\mathbf{R\_{f}} = \frac{1}{2(2\Delta L\_{\rm lf} + L\_{\rm c})\sqrt{\varepsilon\_{\rm rc}}} \mathbf{c\_{l\prime}}\tag{1}$$

where Δ*Lle*, *Le*, ε*re*, and *Cl* indicate variation of the electrical length, the phase length, the effective dielectric constant, and the velocity of the light, respectively. Here, the phase length stands for the geometric dimension of the patch along the direction of antenna's line extension and radiation mode. The effective dielectric constant is rephrased by

$$
\varepsilon\_{\rm cdc} = \frac{\varepsilon\_r - 1}{2\sqrt{(1 + 12H/\mathcal{W}\_c)}} + \frac{\varepsilon\_r + 1}{2}.\tag{2}
$$

where *We*, *H*, and ε*<sup>r</sup>* represent the patch's electrical width, constant thickness, and substrate's dielectric constant, respectively. On the basis of Equation (2), the phase length, Δ*Lle*, is denoted by the dimensions:

$$
\Delta L\_{l\varepsilon} = \frac{(W\_{\varepsilon}/H + 0.264)(\varepsilon\_{rc} + 0.3)H}{2.427(W\_{\varepsilon}/H + 0.813)(\varepsilon\_{cd\varepsilon} - 0.258)}.\tag{3}
$$

**Figure 2.** Micro-strip patch (MSP) antenna. (**a**) Composition of the MSP antenna; (**b**) parts of the patch and the substrate.

The Rf fluctuation depends on the path of the tensile strain (ε*L*). While an observed Rf that is higher than the initial Rf implies the strain was employed in the electrical width direction (*Le*), the value of the Rf should be less under the strain that is applied along electrical length direction (*We*). The frequency oscillation is caused by Poisson's effect by

$$\mathcal{W}\_{\varepsilon} = \left(1 - v\_p \varepsilon\_L\right) \mathcal{W}\_o \quad \text{and} \quad H = (1 - v\_s \varepsilon\_L) H\_{o\prime} \tag{4}$$

where *vp*, *Wo*, *vs*, and *Ho* denote the Poisson's ratios of the patch, the original width, the Poisson's ratios of the substrate, and the thickness, respectively. By a variation of the tensile strain, (*We*/*H*) is elucidated as independent in the case where the proportion of width to thickness is less than one. Here, the Poisson's ratios, *vp*, and *vs*, will be the same, leading to the same ε*re* and ε*r*. Hence, the Rf when a strain, ε, is applied along the vertical direction will be denoted by

$$R\_f = \frac{c\_l}{2\sqrt{\varepsilon\_r}} \frac{1}{(1-\varepsilon)(L\_\varepsilon + 2\Delta L\_{l\varepsilon})} \approx \frac{R\_{o\cdot f}}{1-\varepsilon} \approx R\_{o\cdot f}(1+\varepsilon),\tag{5}$$

where *Ro f* indicates the Rf that is estimated by *Wo* and *Ho*. In this study, an example of the tensile strain applied along the width direction is demonstrated. The patch and substrate are able to derive significant effects on the antenna performance. Existing research has suggested ways to intensify the productivity of the antenna by modifying copious antenna components [35,36], but there exists no rigorous deliberation in terms of the variation of the substrates thickness so far [7]. Such deliberation regarding the variations is essential because most fabrication of the substrate is manual. Also, the interests in printable or wearable antenna applications make the deliberation of non-uniform surfaces valuable. Therefore, it is fairly necessary to assume that the substrate thickness varies; however, these additional concerns demand excessive growth in the number of corresponding responses and DVs. When the engineers decide to weight the additional intricacy, it is typical to encounter the curse of dimensionality, which demands the DR methods in the design process.

Copious DVs result in the curse of dimensionality, which triggers unreliable data assessment in classification and function approximation [37]. Furthermore, the multivariate variables are strongly correlated and arouse the redundancy of the data. The higher redundancy imposes negative impacts on the estimation of the system performance because the inter-correlated data is probably required to represent properties already taken into account. To describe the redundancy, Entropy [38] is often exploited. The uncertainties of a certain random variable are measured by Entropy (*H*) exploiting a probability distribution [39]. Here, *x* referring to a random variable *H* is represented by the following:

$$H\_X = \sum\_{i=1}^m p(\mathbf{x}\_i) \log\_2 \frac{1}{p(\mathbf{x}\_i)},\tag{6}$$

where *p*(*xi*) means the marginal probability of *x*. According to the *H*, the redundancy (*Re*) is able to be re-written as follows:

$$R\_{\varepsilon} = \log\_2 N - \sum\_{i=1}^{n} p(\mathbf{x}\_i) \log\_2 \frac{1}{p(\mathbf{x}\_i)},\tag{7}$$

where log2 *N* indicates the maximum *H* with the entire samples (*N*). In order to elevate the assessment accuracy of the system performance, the DR approach is conducted in this research. The DR technique is generally conducted by the main methods: feature extraction (FE) [40] and feature selection (FS) [41]. While FE transforms the initial features located in the higher dimension to lower-dimensional new features, FS takes the most substantial subgroup of the raw features. Details of the DR methods will be explained in the following sections.

#### *2.2. Feature Extraction (FE)*

FE builds distinct features in the lower dimension using a conversion of the initial features. This DR method can be interpreted as a mapping, which reduces the *Re* of original data. The newly extracted features include most related properties from the initial data. Two major methods, principal component analysis (PCA) [26] and auto-encoder (AE) [42], are explained in this research.

#### 2.2.1. Principal Component Analysis (PCA)

To draw independent and uncorrelated features from heavy correlation, PCA employs orthogonal transformation, which is represented by *Y* = *PX*. The transformed matrix, *Y*, is explained by a transformation matrix, *P*, and the original data set, *X*. In order to conduct a linear transformation, *X* and *Y* will contain numerous examinations and variables, which can be described as *m* and *n*, respectively. Specifically, PCA requires a decomposition of eigenvectors [26] to truncate the data dimension. A covariance matrix, *CY*, can be acquired by

$$\mathbb{C}\_{Y} = \frac{1}{n-1} P A P^T = \frac{1}{n-1} P (P^T D P) P^T = \frac{1}{n-1} D. \tag{8}$$

In Equation (8), *A* = *XXT*, which implies a symmetric condition. *A* consists of the eigenvectors in the rows of the matrix, *P*, and is denoted by *PTDP*. Here, *D* stands for a diagonal matrix that should be connected to equivalent eigenvectors of matrix *A*. The property of an orthonormal matrix allows the inverse of the matrix *P* to be equal to its transpose. Thus, principal components corresponding to *X* will act as the eigenvectors of *PCx* = (*n* <sup>−</sup> 1)−1*XXT*, and the diagonal variables of *CY* should be the *X's* variance.

#### 2.2.2. Auto-Encoder (AE)

While PCA sometimes faces its limitation in non-linear data analysis [26], auto-encoder (AE) is not restricted by the drawback of data linearity. As a special case of the artificial neural network (ANN) [42], AE handles the DR by converting the raw data. AE generally includes three different layers, which are the input, output, and hidden layer, respectively. AE also contains a reconstruction process of the raw data in the output layers. The data will be transformed with a bias and a weight function in the hidden layer. To reduce the gap between the raw and newly recreated data, AE tries to minimize the mean squared error corresponding to the reconstruction [43]. The reconstruction error (*re*) is evaluated by the squared-error cost function,

$$r\_c = 0.5 \| \sum\_{k=1}^{n} \mathcal{W}\_k^T \left( \sum\_{k=1}^{n} \mathcal{W}\_k \mathbf{x}\_k + b \right) + b^T - X\| \tag{9}$$

where *W*, *T*, *x*, and *b* refer to the weights, transpose of the respective vector or matrix, input units, bias, and input vector, respectively. Here, *X* stands for input vector. A sparsity constraint exploits to minimize the total number of hidden units, so that AE can be conducted as the DR. Owing to the constraint, AE has a bottleneck shape. Some of AE's neurons will be inoperative by applying the sparsity constraint to the hidden units. The hidden layer's average activation can be defined by

$$\Psi\_j = \frac{1}{n} \sum\_{i=1}^n \left[ a\_j^{(k)} (\mathbf{x}^{(i)}) \right],\tag{10}$$

where *aj (k)* denotes the activation of hidden layer *k* and the hidden unit *j*. To promote the sensitivity accuracy, an extra penalty term (*P*) is considered.

$$P = \sum\_{j=1}^{s} \rho \log \frac{\rho}{\Psi\_j} + (1 - \rho) \log \frac{1 - \rho}{1 - \Psi\_j} \,\prime \tag{11}$$

where *s* and ρ represent hidden layers and the sparsity parameter, respectively. When the activation value becomes zero, a relation of Ψ = ρ is acquired.

#### *2.3. Feature Selection (FS)*

Unlike FE, FS does not take into account the same dimension when the initial data are reconstructed. It just focuses on grabbing the most significant subgroup in the lower dimension. The importance of the elected subset may not be ensured if the subset contains meaningless properties of the initial data. Diverse studies have been achieved using the mutual information method, genetic algorithm, and single variable classifier method [43], but those methods still have handicaps, that is, huge computational complexity and untrustworthy assessment of the data distributions if the behavior of the distributions is not clearly provided [22]. On the other hand, independent features test (IFT), called the simple hypothesis test [27], is able to quickly discard unnecessary features. For this advantage, IFT will be employed in this proposed framework. The target data are assumed as categorical, so IFT allows all features to be assigned to two categories. This approach is conducted by a scoring value of informative features, SVIF.

$$SVIF(F) = \frac{\left| \mu.(A) - \mu.(B) \right|}{\sqrt{\left(\frac{v.(A)}{n\_1} + \frac{v.(B)}{n\_2}\right)}} > sv$$

From Equation (12), two data sets related to the feature (*F*)'s values can be described as *A* and *B*, where *n1* and *n2* denote the number of the features employed in the categories. With variance, *v*, and mean, μ, the significance value of data can be calculated. In accordance with a threshold, *sv*, indicating

the significance value to eliminate unimportant features, IFT is able to get the best features. In general, it is suggested that the significance value is equal to or greater than 2 [27].

#### **3. DR-integrated Framework for Optimized Design**

On the basis of the developed framework, the DR procedure will be combined with the shape optimization. Even if the framework is expressed with an MSP antenna example, the framework is generalized to be applicable for different design applications, which can be affected by the issue of the curse of the dimensionality. The framework is composed of four steps, as represented in Figure 3. Step 1 shows multivariate input variables that will be randomly created. From Step 2, *Re* of the variables will be minimized by the DR, improving computational productivity. In Step 3, the framework integrated with the DR will find the optimal design. Finally, the reliability of the optimized design of an MSP antenna will be assessed.

**Figure 3.** Proposed optimum design framework with the dimension reduction (DR). ECC, entropy-based correlation coefficient.

#### *3.1. Generation of Multivariate Data Random Thickness*

The design framework in Figure 3 starts with the generation of random multivariate data including a strong correlation in Step 1. For example, in the design space of the antenna application in this study, the correlation will be considered for the substrate thickness. A Gaussian copula is proposed to be used in this study to meticulously create multi-dimensional random behavior. In comparison with the conventional methods that require a linear correlation coefficient [44], the copula takes two clear benefits: (1) it creates multi-dimensional data, although the DVs contain each distinct marginal distribution; and (2) it analyzes the correlation between the random variables that have non-linear behavior. Hence, in Step 1, the most realistic properties in the multi-dimensional space will be demonstrated. A joint distribution, *JXY*, is denoted by the two distinct marginal distributions, *F*<sup>1</sup> and *F*2,

$$J\_{XY}(\mathbf{x}, \mathbf{y}) = \mathbb{C}P(F\_1(\mathbf{x}), F\_2(\mathbf{y})). \tag{13}$$

Here, *CP* represents the copula function (*CP:* [0,1]<sup>2</sup> <sup>→</sup> [0,1]),

$$\begin{split} \text{CP}\_{\rho}(\mathbf{x}, y) &= \frac{1}{2\pi \sqrt{(1 - \rho^2)}} \int\_{-\infty}^{\Phi\_1^{-1}(\mathbf{x})} \int\_{-\infty}^{\infty} \exp\left(-\frac{h^2 - 2\rho hk + k^2}{2(1 - \rho^2)}\right) dh d\mathbf{k}, \\ &\Phi(\mathbf{x}) = \frac{1}{2\pi} \int\_{-\infty}^{\infty} \mathbf{e} \Big(-0.5 \mathbf{s}^2\Big) dh \,\mathrm{\!\_{\prime}} \end{split} \tag{14}$$

where, ρ, *x*, and *y* stand for the linear correlation coefficient, and the marginal distributions of both *x* and *y*. Moreover two different copula parameters are represented by *h* and *k*, and the standard univariate Gaussian distribution is denoted by Φ.

Therefore, the copula function builds the sample dataset by exploiting the parametric multivariate distribution. The samples can then be generated from the copula and can be utilized for stochastic analysis to model and simulate a complex engineering system.

#### *3.2. Entropy-Based Correlation Coe*ffi*cient (ECC)*

Step 2 of the proposed framework involves the DR process. In the high dimensional system analysis, it is required to make an obvious guideline to decide which DR method (i.e., FE or FS) should be accurately employed. In the framework, an explicit criterion is presented to properly select either FE or FS. The criterion can be exploited in accordance with a status of random variables' correlation. As a conventional method, linear correlation estimation [45] measures the degree of the correlation between the data sets. But the approach cannot be supported once features contain non-Gaussian distributions, non-linear correlation, or correlation coefficient of 0.5. To overcome these challenges, it is inevitable to provide a distinct criterion corresponding to an entropy-based correlation coefficient (*ECC*) denoted by *e* [46]. A joint entropy, *HXY*, is calculated by

$$H\_{XY} = \sum\_{i=1}^{n} \sum\_{j=1}^{n} p\_{j\text{int}}(\mathbf{x}\_i) p\_{j\text{int}}(y\_j) \log\_2 \frac{1}{p\_{j\text{int}}(\mathbf{x}\_i) p\_{j\text{int}}(y\_j)} \tag{15}$$

where *pjoint*(*xi*) and *pjoint*(*yj*) are the joint probability distributions of *xi* and *yj*, which are the random variables. As a basic concept of mutual information (MI), entropy can be used to measure dependence, but it does not have a satisfactory scale as the maximum value depends on the size of samples [47]. Hence, the ECC from the normalized MI rescales the MI values to be between 0 and 1. The status of the dependence between the random variables is then calculated by *e*:

$$\varepsilon = \sqrt{\frac{H\_X + H\_Y - H\_{XY}}{H\_X + H\_X}}, \quad 0 \le \varepsilon \le 1. \tag{16}$$

Here, 1 refers to heavy correlations, and vice versa. On the basis of two different ranges of *e*, 0 ≤ *e* < 0.5 and 0.5 ≤ *e* < 1, the proposed framework can provide a guidance to use either FE or FS. The first range indicates that the behavior of the variables is independent and uncorrelated. Thus, the FS will minimize the size of the multi-dimensional data. The second range, on the other hand, means that there is a correlation between multivariate data. Hence, FE will be suitable to reduce the redundancy of the data, *Re*, originated from the strong correlation between multivariate data. When the value of *e* in Equation (16) is close to 0.5, it suggests engineers may employ both DR methods. After the reduction proceeds, the comparison of *Re* values calculated by Equation (7) is required to confirm whether or not sparse features are sufficiently acquired. If the raw features still include higher *Re* compared with the gleaned new features, an additional DR process might be required. Thus, the combination of *e* with *Re* can be considered as a reliable solution that can draw independent and sparse features by reducing the computational cost.

#### *3.3. Optimization*

Once the DR process is conducted in Step 2, an optimization process is then utilized in Step 3 for the optimal design of the application structure or geometry. A function approximation of the actual physical calculations is often required to alleviate the computational burden during the iteration process of the optimization. Therefore, a machine learning-based surrogate modeling technique, namely the artificial neural network (ANN) [27], will be exploited in this research. Compared with the traditional surrogate modeling approaches conducted by linear programming [48], decision trees [49], and discriminant analysis [50], ANN possesses two major benefits: (1) the credibility of ANN is assured in which the decision domain holds complex shapes that are difficult to secure, and (2) ANN is able to be conducted for both classification and function approximation. Essentially, ANN is preferred to be used for the function approximation approach to take the complex contours in the decision domain.

#### *3.4. Predicting Reliability of the Obtained Optimal Design*

Once the DR and the optimization are exploited, the reliability of the acquired antenna's optimum design will be checked. This might be a noncompulsory step based on the ascribed design problem. Generally, when a multi-dimensional system is evaluated, it is frequently computationally exorbitant to analyze further response statistics such as reliability of the system. Hence, although ANN is employed for the function approximation during the optimization process in Step 3, in Step 4, another neural network method, called the probabilistic neural network (PNN), is utilized as a classification process is needed for the reliability assessment. As a major advantage, PNN does not require expensive computational costs because of the fast and straightforward training procedure when compared with ANN. The performance of PNN is assured by the Parzen nonparametric estimator and the Bayes decision rule, which lessen the predicted risk of misclassification [51,52].

The following section presents the applicability of the proposed framework by a design example of the stretchable MSP antenna.

#### **4. Stretchable MSP Antenna-Based Strain Sensor**

In general, antennas that perform in a single frequency band become a barrier of careful surveillance because of the harmful radiation of wireless components and the interruption of data signal transmission. Hence, the efficacy of the single-band antenna cannot to be ensured. For this reason, a dual-band antenna should be utilized to improve the reliability of the stretchable strain-based MSP antenna. In this example, a dual-band antenna functioning at 2.5 GHz and 5 GHz, usually exploited for wireless fidelity, was regarded. On the basis of composition of the MSP antenna, as shown in Figure 2, vital DVs and the material properties are listed in Table 1. Here, μ and *COV* stand for mean and coefficient of variation, respectively.


**Table 1.** Geometric and material properties of the stretchable micro-strip patch (MSP) antenna.

To investigate the behavior of stretchable antennas, a multi-physics analysis that includes both the strain estimation under mechanical deformation of the antenna and the Rf estimation by electrical analysis is required. In the previous research [8], such a multidisciplinary analysis was performed for a stretchable antenna geometry and it was shown that the fluctuating thickness of the substrate resulted in a high return loss. This result was because the uncertainties of the substrate thickness led to the escalation of forged feed radiation and surface waves, which limit the range of the bandwidth. Therefore, in this proposed framework, a design optimization process was also utilized to find the optimal contour of the antenna patch to maximize antenna efficiency by reducing the return loss. The existing geometry in the work of [8] was used as the starting geometry in the optimization process; therefore, first, the dimension reduction and prediction using Steps 1 and 2 of the proposed framework were used for the starting geometry.

In the analysis of the antenna, first, a displacement that results in a tensile strain along the width direction, which is denoted as the *y* direction on Cartesian coordinates, was applied on the antenna substrate, as demonstrated in Figure 4a. In the structural analysis to evaluate the deformation of the antenna under this tensile strain, a linear static finite element analysis (FEA) was employed in this study using the commercial software ANSYS®. As the geometry and the applied displacements are also symmetric, while the material properties are assumed to be isotropic, symmetry boundary conditions were utilized, as shown in Figure 4b. In this study, it was assumed that the maximum applied tensile displacement on the stretchable antenna is 12 mm. The deformation of the stretchable antennas might be larger, but for demonstration purposes, 12 mm was considered in this study as the bound of the displacement of the stretchable antenna. The application can easily be extended for larger displacements in a future study.

**Figure 4.** Schematic of the stretchable MSP antenna. (**a**) Boundary conditions for the tensile test; (**b**) symmetry conditions; (**c**) 202 divided substrate parts containing different thickness and 1 patch with a constant thickness; (**d**) 57 Cartesian coordinates of the patch.

The substrate of the antenna was discretized into 202 rectangular regions that represent the fluctuating substrate thickness, as shown in Figure 4c by white bounds. The outer geometry of the patch shown by gray color was discretized into 57 points, as depicted by white nodes in Figure 4c. Once the deformed shape of the antenna was obtained from the FEA of the antenna model in Figure 4, the return loss and Rf of the deformed antenna geometry were calculated using a commercial software called High Frequency Structure Simulator (HFSS).

#### *4.1. Generation of Multivariate Data for the Varying Thickness*

To model the substrate thickness fluctuation, a stochastic representation of a random field was taken into account. For this purpose, a Gaussian copula was employed to generate the random behavior corresponding to the thickness variation. Figure 5 demonstrates how the copula builds a joint distribution including two distinct sampling sets of thickness variation. First, the roughness of the substrate made of PDMS material was measured by Contour GT-I 3D Optical Microscope. The measured roughness data were added onto the mean thickness of the substrate, which is 0.9705 mm, to describe the variation of the thickness of the whole substrate. In this process, the Gaussian copula generated the random thickness samples engaging several marginal distributions of each sampling set of the thickness variation. Thus, the antenna can be modeled with the substrate, containing 202 small parts that have random different thicknesses and the patch with a thickness of 0.03 mm. Specifically, 500 FEAs were performed with different displacements generated within the range from 0 mm to 12 mm. For each of these simulations, 202 inter-correlated random thickness values were created from the Gaussian copula.

**Figure 5.** The Gaussian copula employed to demonstrate a joint distribution of random sets of thickness variation.

#### *4.2. DR Process for Variation of the Substrate Thickness*

Before conducting the DR, *e* of the initial thickness variation was calculated as 0.684 using Equation (16). This value is greater than 0.5, indicating that the correlation between the thickness data set is strong. This result indicates that FE should be used as the DR according to the proposed framework shown in Figure 3. The performances of two different FE methods, namely PCA and AE, were evaluated in this study. PCA and AE took 32 and 103 dominant components, respectively. In Table 2, the error was calculated for each method by

$$Error = mean \left( \sum \left( \frac{\left| D\_{nuc} - D\_{org} \right|}{\left| D\_{org} \right|} \times 100 \left( \% \right) \right) \right) \tag{17}$$

where *Dnew* and *Dorg* stand for the redundancy of diminished data and that of initial data, respectively. Unlike PCA, AE had higher redundancy because it took extra neuron training processes. However, the dimension taken by AE drew a lower reconstruction error (2.71%) than that taken by PCA (3.27%).

**Redundancy Prediction Coordinates Prediction Error** No *DR* By *PCA* (Redundancy Reduction (%)) By *AE* (Redundancy Reduction (%)) By *ANN* (X coord./Y coord.) By *PCA*/*ANN* (X coord./Y coord.) By *AE*/*ANN* (X coord./Y coord.) *Thickness Coord.* 11.053 5.503 6.127 (44.57%) 7.03%/7.62 % 3.70%/5.41% 4.88%/7.09%

**Table 2.** Redundancy prediction and estimation error of the thickness coordinates. DR, dimension reduction; PCA, principal component analysis; AE, auto-encoder; ANN, artificial neural network.

#### *4.3. ANN-Based Surrogate Model to Predict Antenna Deformation*

(50.21%)

Once the DR was conducted, ANN was constructed to predict Cartesian coordinates of geometry points of the antenna patch, shown in Figure 4d, referring to the shape of the antenna under structural distortion. In the process, the number of variables of thickness variation and simulation was 202 and 500, respectively. The data obtained from these simulations were used to train the ANN model without the DR, with PCA, and with AE. As listed in Table 2, the trained ANN without DR led to a prediction inaccuracy of 7.03% and 7.62% for the *x* and *y* coordinates, respectively. When PCA and AE were used as two different DR methods, on the other hand, this error was reduced, as seen in the results in Table 2.

#### *4.4. Dimension Reduction of Coordinates of the Patch*

In Section 4.2, the DR was employed to analyze the antenna's mechanical behavior to obtain the coordinates of the geometry points of the deformed patch as the outputs. Then, these geometry points are used as the inputs for the electrical analysis of the deformed antenna to calculate the Rf value as the output. On the basis of the inputs, HFSS software was used to check the Rf behavior of the distorted antenna. As output data, 121 resonance frequencies were taken into account.

As this is a new analysis after the mechanical analysis, it was inevitable to conduct the DR for the electrical system behavior as well. The ECC value, *e*, was calculated as 0.442 for the deformed antenna inputs, which are the coordinates of the patch geometry points. This value was less than 0.5, indicating that FS should be used as the DR method in the electrical analysis step. The FS method used in this study was the IFT method and unnecessary coordinates were eliminated based on the significance value of 2. As a result of the IFT method, 28 meaningful coordinates were selected, as shown in Figure 6. Although the *Re* value for the initial coordinates was calculated as 8.382, it was calculated as 5.623 for the deformed coordinates. Thus, it indicates that IFT precisely achieved the *Re* reduction of 32.11% and the correlation of the original coordinates was changed to the independent and uncorrelated ones. In Figure 6, the white points represent the 28 points identified by the IFT method. It is seen that the number of input geometry points for the electrical analysis was reduced from 57 to 28. Although the total number of these selected points is still high, those are the necessary ones to accurately represent the geometry.

**Figure 6.** Selected new significant coordinates.

#### *4.5. Optimization*

Once the analyses were conducted for the initial geometry and the necessary geometry points of the antenna patch were determined, a shape optimization was utilized to improve the performance of the antenna. As explained in the introduction, diverse studies on the strain patch antenna have already been conducted by optimization, but most of them have concentrated on the sizing optimization to take updated height, width, or length of the antenna systems. The optimization may fail to boost the productivity of the antenna because the optimized design almost maintains its original shape. Topology optimization has also been employed to develop a new design, but its efficiency is only considerable under expensive additive manufacturing fabrication. In this study, therefore, a shape optimization was used to find the optimal shape of the antenna patch that minimizes the variation of the Rf and the return loss. It implies that, not only can the possibility of frequency shift induced by structural deformation be decreased, but also the performance of the antenna corresponding to electrical return loss will be maximized. Here, the lower return loss is guaranteed, and the better antenna performance is expected. Unfortunately, however, the explained general equations (Equations (1)–(5)) corresponding to the Rf variation have an assumption that the substrate should have a constant thickness. This assumption only draws a linear relation between the shift of the Rf and the strain. With respect to the stretchable strain MSP antenna including the thickness variation of the substrate, the direction of the Rf shifts could be anticipated by the general equations, but a correct calculation of the shift will not be guaranteed. For those reasons, in this study, a surrogate model to establish a relation between mechanical and electrical behavior of the antenna system was developed, engaging detailed FEA for structural deformation in ANSYS® and electrical analysis in HFSS. The Rf shift and the return loss respecting the structural deformation with additional attention of thickness variation were analyzed with the help of this surrogate model.

To establish objective functions, the substrate thickness variation and distorted patch's shape were considered. The function approximation modeled by ANN was combined with the shape optimization process. To escalate the capability of the antenna, an additional objective function was employed to truncate the return loss.

$$\text{Minimize} \quad \begin{bmatrix} \left| R\_{ff}(\mathbf{x}, f\_{i}) - R\_{fs}(\mathbf{x}, f\_{i\prime}\varepsilon) \right| \\ S(\mathbf{x}, f\_{i\prime}\varepsilon) \end{bmatrix} \quad i = 1, 2. \tag{18}$$

$$\text{Subject to } \left| R\_{ff}(\mathbf{x}\_{\prime} f\_{i}) - R\_{fs}(\mathbf{x}\_{\prime} f\_{i\prime} \boldsymbol{\varepsilon}) \right| \leq \theta\_{\prime} \tag{19}$$

$$S(\mathbf{x}, f\_{i\prime}\varepsilon) + \gamma \le 0,\tag{20}$$

$$\mathbf{x}\_l \le \mathbf{x} \le \mathbf{x}\_{u\nu} \tag{21}$$

$$0 \le \varepsilon \le \Re 0,\tag{22}$$

$$f\_1 = 2.5\,\text{GHz}\,\text{and}\,\,f\_2 = 5\,\text{GHz},\tag{23}$$

where *R*ff and *Rfs* stand for Rf of frequency samples and Rf under applied strain (ε). *x*, θ, *S*, and γ refer to a vector of the DVs, bandwidth regarding a return loss of −10 dB, the return loss, and acceptable return loss, respectively.

The bandwidth corresponding to −10 dB was employed to assure enough antenna operation, and the strain of 30% was established based on the boundary exploited in the FEA. This constraint enables linear static analysis. In order to utilize the benefit of UWB, dual bandwidths, 2.5 GHz and 5 GHz, which are preferred in a wireless fidelity (Wi-Fi) were considered. Moreover, a spline curve was used to model the geometry of the patch antenna with 28 geometry points. The spline coordinates were considered as the DVs, *x*, and the minimum and maximum boundaries of the coordinates were determined to bring optimum shape change of the initial MSP antenna. Particularly, the coordinates of feedline in the length direction were firmly restricted to keep the initial shape of the MSP antenna. The optimization process includes significant attention associated with the strain applied along the width direction. On the basis of the equation of the frequency shift, even though the strain occurs along the length direction, which can make the patch's length varied owing to the Poisson's ratio, the variation will not bring a major impact on the variation of the Rf. Thus, design space around the patch had much more design freedom to improve antenna's operation than that around the feedline. To conduct multidisciplinary design optimization, genetic algorithm (GA) was exploited in accordance with the developed ANN-based surrogate model. For the algorithm, 1000 iterations, crossover probability of 0.96, and mutation probability of 0.01 were used. After conducting the optimization process, a butterfly-shaped antenna was garnered, as shown in Figure 7. Locations of the acquired spline coordinates by the optimization were slightly modified to have a smooth contour. As anticipated, it focused on reducing the variations of the width of the antenna patch, not on that of the length. It could boost the antenna's operation under the applied longitudinal strain. The optimization process drew a reasonable outcome, coinciding with the equations mentioned above in Section 2.

**Figure 7.** Design of the MSP patch antenna. (**a**) Initial antenna design; (**b**) optimal antenna design.

To assure the exactness of the optimized butterfly-shaped design, HFSS was utilized to evaluate the frequency shift. The results of the undeformed antennas are shown in Figure 8. In accordance with the employed bandwidth of −10 dB, an acceptable Rf range of the antenna was assigned from 2.3454 GHz (a) to 2.6539 GHz (b) and from 4.9130 GHz (c) to 5.1070 GHz (d) under 0 mm displacement, as depicted in Figure 8.

**Figure 8.** Rf comparison of initial, deformed, non-deformed optimal antenna, and deformed optimal antenna. (**a**) 2.3454 GHz; (**b**) 2.6539 GHz; (**c**) 4.9130 GHz; (**d**) 5.1070 GHz.

First, the capability of the suggested optimum design was compared with the initial antenna design when the deformation exists for both. Once the deformation of 12 mm is applied to the initial antenna, Rf drastically shifts, ranging from 2.5 GHz to 2.73 GHz and from 5.05 GHz to 5.23 GHz, which states that the antenna functionality will not be ensured because they deviate from each acceptable Rf range regarding 2.5 and 5.0 GHz. The return loss also decreased by 27.6% and 31.8%, respectively. As expected, the performance of the antenna was deteriorated under the deformation. Especially the antenna design in 5 GHz should be restricted owing to miserable return loss, which is higher than −10 dB. Contrary to the initial antenna, on the other hand, the optimized antenna excluding the DR process acquired credits for its flexible possibility under the tensile strain by enabling the Rf to be maintained in the reliable range (θ) and the return loss to be kept under −10 dB, as elucidated in Table 3.


**Table 3.** Comparison of resonance frequency (GHz) and return loss (dB) under strain test.

Also, a comparison of two additional designs regarding the DR process antenna was conducted. In Table 3, the design of the initial MSP antenna design including DR led to a better operation than that of the initial antenna itself. Although, its improvement is not enough to dominate the initial antenna design with an optimization. As predicted, the optimized antenna design with DR had notable improvement when compared with all other cases. Therefore, the uncertainty-reduced optimal MSP antenna design could be a dominant approach because general wireless receivers would be employed without any additional extension work on the bandwidth.

#### *4.6. PNN for Classification of the Resonance Frequency*

In the final step, the reliability of the optimal MSP antenna was analyzed in accordance with the reliable range of θ. The antenna must operate within the range even though it is under deformation. The range was estimated by the non-deformed MSP antenna. From the previous step, the dual bands antenna under no displacement took the following range (θ):

$$\begin{aligned} \text{2.3454 GHz} & \le \theta\_{f\\_2.5\text{GHz}} \le \text{2.6539 GHz} \\ \text{4.9130 GHz} & \le \theta\_{f\\_5.0\text{GHz}} \le \text{5.1070 GHz} \end{aligned} \tag{24}$$

For reliable classification, the limit state function was exploited. According to the concept of the limit state function [53], the capacity or resistance of the antenna must be appointed to 2.5 GHz. However, if the absolute value of *g* is greater than 0.3085 (difference between 2.3454 GHz and 2.6539 GHz) and 0.194 (difference between 4.9130 GHz and 5.1070 GHz) at 2.5 GHz and 5 GHz, respectively, the designed strain MSP antenna will stay in class *B*, implying that the system must be ignored owing to the unreliable Rf. The Monte Carlo simulation (MCS) [53,54] was employed to estimate the probability of failure (*Pf*). On the basis of the *Pf*, the results of MCS with 10,000 random variables and PNN with 121 random variables were compared to assure results of the reliability estimation. As explained in Table 4, each *Pf* difference between PNN and MCS was 6.46% and 6.91% at 2.5 GHz and 5 GHz, respectively. In comparison with the original variables, the truncated coordinates obtained by IFT had a *Pf* of 0.291 and 0.325 by drawing 4.67% and 4.06% of *Pf* increase, respectively. Furthermore, the *Pf* acquired by IFT was close to that obtained by MCS (9.48% and 9.41%, respectively), which is less than 10%, so it can be a reliable threshold, eliminating computation-intensive tasks, for the classification process.


**Table 4.** Probability of failure conducted by probabilistic neural network (PNN) and Monte Carlo simulation (MCS).

#### **5. Conclusions**

In this research, a design framework exploiting a DR approach, machine learning-based surrogate modeling, structural optimization, and reliability assessment was proposed to handle a multivariate and multidisciplinary engineering system. The efficiency of the proposed design framework was addressed by a stretchable MSP antenna-based strain sensor.

Compared with the conventional rigid antenna, the efficiency of the proposed stretchable antenna design was highlighted because it contains careful consideration for all possible mechanical flexibility, which weights realistic system prediction and estimation. Such consideration demands meticulous examination on the flexible substrate that includes non-uniform thickness. In this research, therefore, the non-uniformity was regarded as the DVs to represent geometric uncertainty captured by a Gaussian copula function. With the copula function, a non-uniform substrate thickness was represented by 202 subparts, taking the mean of 0.9705 mm and the COV of 0.05. However, a huge number of DVs escalate the redundancy and complexity of data, inhibiting precise system prediction and estimation. In order to resolve the issue, the proposed framework employed a DR process, particularly FE and FS. In the DR process, ECC (*e*) was exploited as a clear guideline that suggests a better DR method grounded on data behavior. On the basis of the ECC estimation, FE was employed to reduce data redundancy regarding

the substrate non-uniformity in which *e* is greater than 0.5, whereas FS eliminated unnecessary coordinates of the patch within the specific case (0.5 ≤ *e* < 1). Both FE (PCA and AE) and FS (IFT) derived good redundancy reduction of 50.21%, 44.57%, and 32.11%, respectively. With the criterion, engineers can make a decision on the applicable DR process and management of multivariate data.

Moreover, owing to the limitation to formulate a relationship between structural deformation and electrical antenna's response, the ANN-based surrogate model employing the multivariate data purified by the DR process was developed to predict a complex engineering system. As the second drawback of the stretchable MSP strain sensor, the narrow bandwidth that restricts the functionality of the sensor was handled by a structural shape optimization. Compared with the conventional antenna shape, the proposed optimum shape drew the antenna's performance improvement of 5.77% at 2.5 GHz and 29.03% at 5 GHz within the reliable frequency range. A new optimum design implementing the DR process also escalated the performance improvement of 29.05% at 2.5 GHz and 35.08% at 5 GHz. Thus, it is confirmed that the developed optimal design maximized the frequency stability within a reliable bandwidth and antenna performance under structural deformation. In the final step, the reliability of the stretchable antenna for the strain sensor was assessed by PNN and validated by MCS. PNN evaluated the efficiency of the proposed design by showing the classification accuracy of 96%. The accuracy of PNN was validated by that of MCS. The results obtained for the stretchable strain MSP antenna show that the proposed design framework with the uncertainty characterization and dimension reduction is effective on multi-physics-based and multi-objective design processes.

**Author Contributions:** Conceptualization, S.-K.C.; methodology, S.H.; software, S.H.; validation, S.-K.C. and R.M.G.; formal analysis, S.H.; investigation, S.H.; resources, S.-K.C.; data curation, S.H.; writing—original draft preparation, S.H.; writing—review and editing, S.-K.C. and R.M.G.; visualization, S.H.; supervision, S.-K.C.; project administration, S.-K.C.; funding acquisition, S.-K.C. and H.-J.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially funded by Chung-Ang University Research Grants in 2016.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Product Service System Availability Improvement through Field Repair Kit Optimization: A Case Study** †

#### **Eun Suk Suh**

Graduate School of Engineering Practice, Institute of Engineering Research, Seoul National University, Seoul 08826, Korea; essuh@snu.ac.kr; Tel.: +82-2-880-7175

† This article is a re-written and extended version of "Improving Product Service System Availability through Field Repair Kit Optimization" presented at the 2nd East Asia Industrial Engineering Workshop (EAWIE 2015), Seoul, Korea, on 7 November 2015.

Received: 11 September 2019; Accepted: 9 October 2019; Published: 12 October 2019

**Abstract:** Product service system (PSS) is becoming a popular business model, where companies offer product based service to customers to realize steady recurring revenue. However, to provide PSS-based service to customers in reliable way, PSS need to be supplemented with a field repair kit onsite, in case of parts failure and PSS shutdown. The field repair kit consists of frequently used spare parts in multiple quantities. However, mismatch in spare parts type and quantities in the field repair kit will results in sub-par performance of PSS for both customer and company. In this paper, a case study involving industrial PSS repair kit optimization is presented. In the case study, the field repair kit for complex industrial printing system is cost optimized, while satisfying the system availability requirement, specified by the maintenance contract between the company and the customer. Key analysis steps and results are presented to offer insight into the PSS field repair kit optimization, offering useful references to industrial practitioners.

**Keywords:** product service system (PSS); availability; field repair kit

#### **1. Introduction**

Research on the product service system (PSS) [1,2] is increasing as many leading global industries are shifting their focus toward product based services. For example, in the aircraft engine industry, the old business model aimed to sell an engine to an aircraft manufacturer and provide engine maintenance, if necessary. Over time, this business model has evolved into a service centric model, often called the "power by the hour" approach [3], where aircraft engine companies own the engines installed on customer's aircrafts, and charge their customers for actual flight hours. This paradigm shift impacted many facets of product development and lifecycle management, including product requirement definition, subsystem design, service development, and total product lifecycle cost analysis.

According to general literature survey by Beuren [1], academia defines PSS as "a combination of products and services in a system that provides functionality for consumers and reduces environmental impact" [4]. Similarly, Baines et al. [5] defined PSS as "an integrated product and service offering that delivers value in use to the customer". Tukker [6] further categorized PSS into eight different types, which are clustered into product-related services, use-oriented services, and result-oriented services. The PSS presented in our case study is the product-related service type, whose company offers a product and related services (e.g., financing, maintenance contract, spare and consumable parts supply) throughout the lifecycle of the product.

There are four different essential elements in the PSS: product, process, related human role, and service. Product typically consists of actual hardware equipment used to provide specific functions

required. Process is the sequence of actions required to provide the desired service. A related human role consists of activities that need to be performed by personnel who are part of the PSS. Finally, the service is the output of all these equipment and activities that is provided to customers. One such example of PSS is an amusement park ride. The product is the park ride itself, which is necessary for providing necessary entertainment service to ride users. The process is a sequence of activities necessary to provide ride service, which consists of getting customers into the ride, operating the ride equipment, getting customers out of the ride, and equipment maintenance. Related human role include operating and maintaining ride equipment. The service aspect is the ride experience itself, enjoyed by customers.

For PSS, one of the key performance metrics is the system availability. The PSS must be operational to meet the availability requirement promised to customers. For the amusement park ride example previously mentioned, the ride apparatus must be operational during park business hours to provide customers with a satisfactory ride experience. Not meeting this availability requirement, due to ride apparatus failure, will result in loss of customers and revenue. One way to maximize the uptime availability of the PSS is to keep a field repair kit, consisting of spare parts for the PSS, on the customer site, in case of such a failure. In many cases, product design teams rely on their prior experiences to determine the types and quantities of spare parts in field repair kits. However, this often result in a sub-optimal field repair kit, containing some spare parts that are never used, or carrying less-than-necessary quantities of some spare parts that are always used. This mismatch of spare parts inventory may cause unacceptable downtime and lost revenue for PSS customers and providing companies.

In this paper, a PSS field repair kit optimization case study is presented in detail. In the case study, a field repair kit for a complex PSS (industrial printing system) was optimized in terms of field inventory kit cost, while satisfying the availability requirement set by contract with the customer. A high fidelity simulation PSS simulation model was created to simulate spare parts usage during the PSS operation. Using the model, a cost optimal field repair kit was identified. Subsequent analysis was performed to determine the confidence level that the PSS is capable of achieving the imposed PSS availability requirement. The paper is organized as follows. A survey of related research literature is presented in Section 2. The case study overview, optimization process, analysis of the results, and discussion are presented in Section 3. The paper closes with conclusions in Section 4.

#### **2. Literature Review**

Traditionally, product manufacturers have focused their primary efforts on product development and sales. For these traditional manufacturing firms, the term "service" consisted of maintenance and repair of their products in the field as needed. However, as companies seek ways to create recurring revenues throughout the lifecycle of their products, the trend to integrate the product with the services offered to the customer has increased. This approach has advantage of continuous and stable revenue generation for the company because it establishes a long term relationship between the company and the customer [7], as well as sustaining environment from the social and company perspectives [8,9]. There have been successful examples of such practices, which include Gage Products and PPG Industries [10], where companies have integrated their core products with the total service offered to the customer. For example, major printing systems manufacturing companies, such as Xerox and HP, offers total document management services to various customers to optimize their total document production workflows.

Once the PSS is deployed in the field, it is of critical importance that the deployed PSS is operational, meeting the contracted availability requirement promised to the customer. One way to meet the availability requirement is through implementing appropriate maintenance policy, which includes preventive maintenance and predictive maintenance. There have been several works published in the field of preventive maintenance and predictive maintenance. Relevant works on preventive maintenance include work by Adhikary et al. [11], who proposed multi-objective genetic algorithm

approach to optimize availability and maintenance cost through preventive maintenance scheduling model, work by Moghaddam [12], who proposed a nonlinear mixed-integer optimization model for a manufacturing system, and the work by Mokhtar et al. [13], who proposed a maintenance policy optimization framework using Bayesian networks. Additionally, related works on predictive maintenance include work by Van Horenbeek and Pintelon [14], who proposed a dynamic predictive maintenance policy and compared it with other various maintenance policies to show that the proposed policy results in significant cost reduction. Another work by Wang et al. [15] proposed cloud-based predictive maintenance paradigm to advance the state of intelligent manufacturing. Finally, there is a research trend that aims to link big data to improve predictive maintenance [16,17].

Another way to achieve required PSS availability is through implementing component redundancy and keeping an onsite spare parts field repair kit, consists of parts that are expected to fail [18]. This will prevent the unintended PSS downtime due to lack of critical spare parts. However, deciding on the quantity of each spare part of the PSS field repair kit requires information on the total PSS life cycle, average usage per designated period, and individual spare part's reliability. The management of spare parts inventory is a well-developed research topic. Published works have investigated topics such as the allocation of spare parts inventory within a multi-echelon supply chain [19–21], spare parts inventory and reliability decision framework with service constraints [22,23], combined optimization of preventive maintenance and spare parts inventory [24,25], and obsolescence management [26]. Recent advances in PSS-related research has expanded to other important areas, such as sustainable product service system design [27,28], cloud based product service system design [29,30], product service system implementation for the smart city [31], and incorporation of digital twin concept to the product service systems [32]. Recent works on field repair kit optimization include the joint planning and optimization of spare parts inventory and service engineer staffing [33,34].

In this paper, the featured case study focuses on the cost optimization of the onsite spare parts repair kit, subject to the PSS availability requirement specified by the contract between the company and the customer. Once the field repair kit is optimized for a specific level of availability requirement, additional analysis will determine the level of confidence of the PSS in meeting the availability requirement with the repair kit. The work presented in this study contributes to existing literature by introducing a real industry PSS case study that can serve as a reference for other industry practitioners.

#### **3. Product Service System Case Study**

#### *3.1. Case Study Overview*

Xerox, one of the world's leading printing system manufacturers, has a fleet of industrial printing systems deployed in the field. Industrial printing systems are installed at customer's print shops, where they become part of production systems to produce printing goods, which, in turn, are sold to their customers. Xerox (the company) and the customer typically have a recurring fee based maintenance contracts, in which the company maintains customer's printing system on regular basis, supplies necessary consumable parts required for operation, and performs repairs in case of printing system failure. Usually, as part of the maintenance contract, the company is responsible to have its printing system available for production above a certain agreed threshold during the regular operating hours. Failure by the company to meet the contracted system availability requirement may result in penalty to the company, or even worse, cancellation of the maintenance contract by the customer. In that sense, the industrial printing system can be classified as a PSS for Xerox, and will be referred to as such in subsequent sections.

Since unexpected system failures are inevitable, the company supplies customers with a field repair kit, consists of frequently replaced spare parts. Typically, the composition of the spare parts and the quantity of each spare part in the field repair kit are determined by printing system design team, based on historical performance of similar parts in other previously launched products. After the

product launch, the composition and quantities of spare parts in the field repair kit are updated as the true life of spare parts become known.

However, on many occasions, this result in a sub-optimal field repair kit, where some parts are never used due to less than expected parts failure, while some parts are in constant shortage and demand as they are used more frequently than expected. This mismatch creates an undesirable situation for both the company and customers. Customers experience more unexpected PSS downtimes, above the threshold established by the maintenance contract, due to a shortage of critical spare parts. This will impact a customer's ability to meet their business goals. As for the company, mismatch in field repair kit spare parts type and quantity may result in their failure to meet the availability requirement promised to the customer, jeopardizing their credibility, and sometimes terminating their contract with that particular customer. Another issue concerns the field repair kit inventory itself. In some cases, the company financially owns the field repair kit inventory, tying up corporate capital in the form of spare parts inventory. If some spare parts become obsolete while in a company's possession, the company would incur unnecessary financial loss due to the extra cash that was tied up in the form of unused parts, and the additional cost for management of obsolete parts. The constant shortage of some spare parts can create an extra cost issue in terms of expedited shipment from parts suppliers and to customer's sites. In this case study, a simulation based optimization framework was created and utilized to configure a cost-optimal PSS field repair kit to satisfy the specific availability requirements. One of several industrial printing system produced by the company was used as a specific example.

To optimize the PSS field repair kit and to perform further analysis, following steps were followed. First, a high fidelity PSS simulation model was created to simulate the system operation and availability, subject to spare parts failure. Next, using the simulation model, a series of key parameter variation experiments were conducted to yield PSS availability under specific parameter settings. The obtained results were used to construct a regression based PSS availability model that would describe PSS availability as the function of key parameters. Regression model was then used to optimize the field repair kit in terms of the inventory cost, while satisfying the availability requirement. Finally, a Monte Carlo simulation of PSS with the optimized field repair kit was performed to estimate the confidence level for the PSS to meet the specific availability requirement.

#### *3.2. Case Study Assumptions*

Several assumptions were made to establish the boundary within which the case study is valid, and to reflect reality. Simulations, subsequent optimizations, and analysis were performed within the framework of the assumptions.

• Definition of availability: in this case study, the availability of the PSS is defined as follows. When the PSS shuts down due to a spare part failure, the part is replaced with a fresh part from the field repair kit. If there is no replacement part present in the field repair kit inventory, it must be shipped from the central warehouse. This is recorded as a PSS down incident. The availability metric is the percentage of all PSS part requests that is fulfilled by immediate replenishment from the field repair kit. For example, if there were 100 requests for various spare parts, and 95 of them were fulfilled from the field repair kit without placing emergency order from the central warehouse, the availability of the PSS is 95%.

• PSS and production system configuration: for the case study, the production system consists of two PSSs of the same type in a serial configuration. The field repair kit optimized in the case study is configured to service both PSSs.

• Field repair kit spare parts composition and quantities: the design team provided the list of spare parts in the field repair kit, along with their respective quantities and the estimated life. For the case study, all parts identified by the design team were included in the optimization. A total of 55 spare parts were included in the field repair kit.

• Periodic PSS usage: the PSS usage is defined as the number of prints produced per time period (month in this case study). For the case study, the nominal usage is set to 8.0 million, as provided by the design team responsible. The usage data was derived from past historical data of customers who have used similar type of PSS. Additionally, it was assumed that all spare parts fail based on PSS usage, not based on time.

• Onsite conditions: availability of PSS also depends on other aspects, such as experience of PSS operation and maintenance staff and production equilibrium. It was assumed that the level of expertise of operation and maintenance staff was competent to upkeep the PSS at the desired performance level, as it typically is for real onsite facilities. Additionally, it was assumed that the production equilibrium is maintained.

• Simulation assumption: for the Monte Carlo simulation, spare parts life distribution and periodic PSS usage distribution were provided by the design team. Additionally, simulation time was set to be equal to the typical contract term agreed between the company and the customer. For this case study, the contract term was set to five years, with a single simulation period equal to one month.

#### *3.3. Product Service System Simulation Model*

The first step in the proposed process is to create a high fidelity PSS model. The Anylogic™, software was used to create the PSS model to simulate PSS usage and usage-based spare parts failures. Figure 1 shows the basic structure of the PSS discrete event simulation model. The PSS selected for the case study is typically sold in pairs (quantities of two) to customers, and used in serial configuration to process customer orders for the duration of the contract period. There are work in process (WIP) product inventory as shown in the model. When the simulation starts, document orders start to arrive into the facility. Each order contains a request for certain quantities of a finished product, which, in this case, is the document produced in the desired quantity. The production line initiates product manufacturing, and utilizes the PSS that is part of the total production system. Once requested product quantity is manufactured, the order is complete and the finished document products are processed for delivery to final recipients. The simulation model records PSS usage, which is equal to the total quantity of the products manufactured. When the PSS usage reaches a threshold for a particular spare part with a specific usage based failure rate, the part fails, and the PSS shuts down. The failed part is replaced with the part from the field repair kit inventory. Once the spare part is replaced, the usage count for that specific spare part is reset and the production resumes. The replaced part is replenished from the central warehouse to fill the field repair kit. If the failed part is out of stock at the customer's site, then an emergency order is placed, and the part is shipped from the central warehouse immediately, using one-day express shipment. When this happens, it is recorded as a down incident for the PSS.

**Figure 1.** Simulation model flow chart for industrial printing product service system.

Figure 2 shows the Input-Process-Output (IPO) diagram for the PSS simulation model shown in Figure 1. Required input data are spare part life, quantity of spare part onsite and the monthly usage of PSS. Once the simulation ends, the model yields three specific outputs: (1) total demand for the spare part required by PSS, (2) the total number of PSS down incidents due to spare part shortage, and (3) the PSS availability.

**Figure 2.** Input-Process-Output (IPO) diagram of the simulation model.

#### *3.4. PSS Availability Model*

The second step in the overall process is to construct the PSS availability regression model by performing key parameter variation experiments for PSS using the simulation model. A model of PSS availability as a function of input parameters shown in Figure 2 is created. Table 1 shows the key parameters and ranges of the parameter values used for the case study.

To create a reliable model, 240 different experiment runs were set up and performed using different deterministic parameter values, and resulting PSS availability values were recorded. Table 2 shows results of the 15 deterministic experiment runs out of 240 runs performed. For example, experiment number 2 in Table 2 shows a spare part with a life of 10 million uses (0.1 failure per million use), with only one provided in the field repair kit. When the simulation was performed with PSS monthly usage of 10.2 million, the availability of the PSS overall was 15% due to shortage of the part. Availability for each run was recorded. Once all experiment runs are complete, a regression based PSS availability model was constructed using parameters shown in Table 1. The regression model provides a good surrogate means to calculate PSS availability using spare part failure rate, onsite spare part quantity, and PSS usage per period. One thing to note is that the regression model is valid within parameter value ranges shown in Table 1.


**Table 1.** Key parameters and their value ranges for simulation and optimization.

**Table 2.** Partial results of deterministic PSS simulation.


#### *3.5. Cost Optimization Model*

For the next step of the process, the PSS availability regression model created from the previous step is embedded into the field repair kit cost optimization model. The IPO diagram for the optimization model is shown in Figure 3.

**Figure 3.** IPO diagram of the field repair kit cost optimization model.

For the field repair kit cost optimization, several input parameters are required. The first parameter is the required PSS availability target. This is typically dictated by the customer who utilizes the PSS when the maintenance contract is signed. The second parameter is the unit cost for each spare part in the field repair kit. The third parameter is the quantity of each spare part in the field repair kit, and the last parameter is the failure rate of each spare part. Table 3 shows a partial list of the actual spare parts included in the field repair kit, with actual parameter values. It should be noted that the unit cost of each spare part is added, since the objective of the optimization is to minimize the total cost of the field repair kit. For the monthly PSS usage, it was set at 8.0 million, based on field data for similar systems. Optimization was performed using spreadsheet-based commercial software (QuantumXL™), which uses a heuristic algorithm.


**Table 3.** Partial list of spare parts and quantities in the PSS field repair kit originally proposed by the design team.

For the optimization, three different scenarios were formulated after discussion with the design team. Table 4 lists proposed scenarios. The first scenario serves as the base scenario, where the PSS operates with the original field repair kit composition shown in Table 3. In the second scenario, the PSS availability requirement parameter is set to 95%, and the field repair kit is optimized for total cost. For the third scenario, the PSS availability requirement parameter is relaxed to 90%.


**Table 4.** Key parameters and their value ranges for simulation and optimization.

#### *3.6. Optimization Results and Discussion*

Table 5 shows a partial list of the optimized field repair kit spare parts shown in Table 3. The list for scenario 1 shows the original quantities of spare parts proposed by the design team. Second column shows the optimized spare parts list for 95% availability requirement. The last column shows a partial list of spare parts, optimized for 90% availability requirement. Observation of the list reveals some key insights. In the original repair kit composition, two units of spare part A were kept on the customer site. However, the optimum composition eliminated spare part A altogether from the field repair kit, mainly because of the part unit cost, which was significantly higher than the cost of any other part in the field repair kit. On the other hand, part J, which has a very low part unit cost, saw an increase in quantity from the original quantity proposed by the design team. The observation of the optimization results is that in balancing the field repair kit cost and the availability requirement, the optimized field repair kit is composed in a way that the PSS will only shut down for a part that is too expensive to be held onsite, and not for a low cost part.


**Table 5.** Partial list of spare parts and their quantities in the optimized field repair kit for each scenario.

Figure 4 shows, for each scenario, the total of the cost of the PSS field repair kit. The results show that in order to meet the PSS availability requirement of 95%, the total cost of the field repair kit must increase by more than 60% over that of the originally proposed field repair kit composition in Scenario 1. However, when the PSS availability requirement constraint is relaxed to 90%, the cost of the field repair kit was significantly reduced, nearly down to 50% of the field repair kit cost shown in Scenario 2. In other words, increasing the availability requirement by 5% resulted in a 50% increase in the total cost of the field repair kit. Additionally, it is observed that if 90% PSS availability is acceptable to the customer, the optimized field repair kit can be composed for lesser cost than the original field repair kit proposed by the design team.

**Figure 4.** Field repair kit cost optimization results for scenarios analyzed.

As in the last step, using the optimized field repair kit composition, Monte Carlo simulation was performed to estimate the confidence level for which the PSS can meet the availability requirement for each scenario. Two parameters were used. The first parameter was the failure rate of each spare part. Distribution of spare parts failure rates were set based on the historical reliability data of the part or similar parts. The second parameter was the PSS usage. Since PSSs are used by many customers, but with different monthly usage, the PSS monthly usage was set with triangular distribution with a minimum value of 5.7 million uses and a maximum value of 10.2 million uses, with a mean usage of 8.0 million per month. This was based on the estimate from the design team. The resulting distribution of availability requirement was analyzed to determine the confidence level for the optimized field repair kit to meet the target PSS availability. The confidence level of availability that the PSS can achieve with the optimized field repair kit is plotted in Figure 5. For each scenario, 10,000 simulation runs were performed with the spare parts failure rate distribution and the periodic PSS usage distribution. The horizontal axis represents the confidence level, while the vertical axis represents the PSS availability.

**Figure 5.** Availability confidence level of the PSS with optimized field kit in place.

In the figure, the topmost curve shows the confidence level of the PSS with the field repair kit optimized for 95% availability. The specific data point singled out shows that this particular PSS with the optimized field repair kit can meet the 95% availability requirement with a 75% confidence level. The second curve shows the confidence level of the PSS with field repair kit optimized for 90% availability. It shows a 70% confidence level for meeting the 90% availability target. Finally, the two optimized field repair kits, in contrast with the original field repair kit in the base scenario, performed particularly well in terms of expected confidence level to meet the PSS availability requirements. Results obtained and shown in Figures 4 and 5 and Table 5 can be a valuable guideline for assessing the PSS availability with a specific field repair kit. Cost savings realized by the adjustment of the field repair kit inventory should be balanced against the expected PSS confidence level for meeting availability requirement promised to customers. The company can also create a mitigation plan based on the confidence level to minimize the risk of not meeting the availability requirement.

#### **4. Conclusions**

In this paper, an industrial PSS field repair kit optimization case study of was presented. The objective of the case study was to cost optimize the composition of a field repair kit placed at a customer's site, while maintaining a specific level of PSS availability. A high fidelity PSS simulation model was created and utilized to construct the PSS availability regression model. The availability model was used with other spare parts information to cost optimize the field repair kit, subject to the PSS availability requirements. Monte Carlo simulation was performed to assess the confidence level at which the optimized PSS field repair kit can meet the availability requirements.

Results showed that the PSS field repair kit can be cost optimized to satisfy the availability requirements imposed by the contract with the customer. It was also revealed that the marginal cost of availability improvement was very high, resulting in an almost twofold increase in cost to improve the availability by 5%. The investigation of the optimized PSS repair kit composition showed that the PSS should never shut down for a low cost spare part, but should only shut down for a very expensive spare part that is too expensive to be included in the field repair kit. Finally, Monte Carlo simulation results showed that the optimized PSS field repair kit performed better than the originally proposed field repair kit, meeting the availability requirements at a higher level of confidence. The overall case study results were encouraging in that the modeling and optimization process shown in this paper is applicable to other complex systems.

Future works regarding the improvement of PSS can expand into multiple directions. The latest advances in big data, information and communications technology (ICT), and internet of things (IoT) can greatly improve preventive and predictive maintenance practice for PSS. This can be accomplished through better prediction of parts failure interval and failure modes through analysis of historical data and real-time feedback from the PSS. Big data and IoT can be utilized to design the improved spare parts supply chain as well, providing information for an appropriate inventory level, stocking locations, and reorder points based on PSS location and usage. Another direction the research can take is the improvement of the PSS design: how can PSS be architected for optimum preventive or predictive maintenance policy? How does the PSS architecture impacts spare parts logistics or its supply chain structure? In addition to the topics mentioned, another fundamental direction that the PSS related research can take is to explore ways to reduce the limitation imposed for implementing PSS, such as willingness to adopt the PSS by companies, willingness to accept the PSS by customers, and dealing with environmental implications. All these are very promising future research topics to explore.

**Funding:** This research received no external funding.

**Acknowledgments:** The work presented in this paper in its entirety was conducted while employed at Xerox Corporation.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


© 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Research on Optimized Product Image Design Integrated Decision System Based on Kansei Engineering**

**Lei Xue 1, Xiao Yi <sup>1</sup> and Ye Zhang 2,\***


Received: 18 December 2019; Accepted: 7 February 2020; Published: 11 February 2020

**Abstract:** In order to facilitate the development of product image design, the research proposes the optimized product image design integrated decision system based on Kansei Engineering experiment. The system consists of two sub-models, namely product image design qualitative decision model and quantitative decision model. Firstly, using the product image design qualitative decision model, the influential design elements for the product image are identified based on Quantification Theory Type I. Secondly, the quantitative decision model is utilized to predict the product total image. Grey Relation Analysis (GRA)–Fuzzy logic sub-models of influential design elements are built up separately. After that, utility optimization model is applied to obtain the multi-objective product image. Finally, the product image design integrated decision system is completed to optimize the product image design in the process of product design. A case study of train seat design is given to demonstrate the analysis results. The train seat image design integrated decision system is constructed to determine the product image. This shows the proposed system is effective and for predicting and evaluating the product image. The results provide meaningful improvement for product image design decision.

**Keywords:** product image design; Kansei Engineering; integrated decision system; qualitative decision model; quantitative decision model; train seats

#### **1. Introduction**

The development trend of industrial design is toward more efficient, intelligent and systematic. Meanwhile, product design pays more attention to the improvement of the user experience. The development of industrial design decision system not only promotes the economic and cultural exchange between regions, but also has a profound impact on the long-term development of various related industries [1]. With the increasing demand for user experience design, user perception centered design will become an important factor affecting product design [2]. Therefore, based on meeting the product function, we should comprehensively improve the user physiological and psychological multi-dimensional emotional needs, optimize the overall product design level, promote the user experience, and raise the product design quality.

Kansei Engineering is a kind of theory based on the design science, psychology, cognition and other relevant disciplines, which can lead human perceptual analysis to the field of engineering technology [3]. Kansei Engineering not only helps designers understand users' perceptual needs, but also optimizes product design process and reduce product design cost. Kansei Engineering is a method of using engineering technology to explore the relationship between the sensibility of user and the design characteristics of product. It can transform the perceptual needs by users into the design elements of products. In the field of product design, Kansei Engineering is a user-oriented product

research and development technology. The modern tools and technologies are used to quantify the user's perceptual information into the design parameters of product, as shown in Figure 1.

**Figure 1.** The Kansei Engineering of product design process.

Image is a kind of the psychological symbol that the products give to the users, and is also the psychological image or concept obtained by association and imagination. The users' cognition for the product, through the perception of five senses, causes users' resonance. Finally, a comprehensive psychological image is formed. Image has certain subjectivity and fuzziness. The product appearance design should not only meet the users' basic aesthetic needs, but also meet the users' perception needs. The basic design elements of product constitute the product appearance design. The design elements of product appearance design generally contain form, color, material, texture, pattern and so on. The formation of product image comes from the users' cognition. The design elements of product are taken as the language of communication with the user. Therefore, the formation process of product image is shown in Figure 2.

**Figure 2.** The formation process of product image.

Wei proposed a hybrid adaptive ant colony algorithm to realize product family multi-objective optimization design through scale-based product platform theory model [4]. Lei presented a Decision Support System (DSS) for market-driven product positioning and design, based on market data and design parameters. The proposed DSS determines market segments for new products using Principal Component Analysis, K-means classification [5]. Yang presented a decision support system based on a bi-level fuzzy computing approach that incorporates qualitative and quantitative product attributes in determining the manufacturability of a product [6]. A product design elements evaluation model was proposed, which was constructed by eye-tracking experiments, using eye movement indicators such as first gaze time, gaze order, and number of times of return, to accurately and effectively analyze product model and quantify users' emotion. The purpose of sorting the product model contribution was achieved through the weight calculation of the design elements [7]. Understanding the affective needs of customers is crucial to the success of product design. Hybrid Kansei Engineering system (HKES) is an expert system capable of generating products in accordance with the affective responses. HKES consists of two subsystems: forward Kansei Engineering system and backward Kansei Engineering system [8]. Wu proposed a preference-based evaluation-fuzzy-quantification method to determine the priority of the development of attractive factors of the product [9].

The state of the art has the following deficiencies: (1) The previous research lacks the product image design qualitative decision by the quantitative research methods; (2) It lacks the integrated decision system of product image design for the purpose of the total image design for assistance of the collaborative work of user, designer and expert; (3) The existing product image design model lacks the discussion of the assessment and verification with the computer-aided and 3-D printing of product design.

In our paper, we extend the research objectives to the product material design and texture design. This research proposes a new method for combining the product image design qualitative decision model and quantitative decision model. With the two models, the optimized product image design integrated decision system is constructed based on the Kansei Engineering experiment. In the process of the performance evaluation, the technology of the computer-aided and 3-D printing assist the evaluation with Kansei image of product image design integrated decision system.

Therefore, this paper explores the product image design system under the research framework of Kansei Engineering. The qualitative and quantitative research on the correlation between product design elements and user Kansei image is conducted, and the qualitative and quantitative decision models of product image design are constructed. Then, the integrated decision system of product image design is improved. The paper introduces the methodology that is divided into four steps, shown in Figure 3, including the establishment of qualitative decision model of product image design, quantitative decision model of product image design, product image design integrated decision system and its performance evaluation. Meanwhile, this paper takes the high-speed train seat design as the research object. Kansei Engineering experiment and design theory research are carried out deeply, and the more scientific, efficient and intelligent product image design evaluation system is put forward.

**Figure 3.** The four steps of the methodology process.

#### **2. Kansei Engineering Experiment**

The Kansei Engineering experiment process of product image design integrated decision system mainly includes design investigation stage and data statistics stage as shown in Figure 4. The design investigation stage is divided into two steps. The first step is to select the experimental subjects, experimental samples, product form samples, color samples, material samples and texture samples. Then the second step is to analyze the design elements of experimental design samples separately, such as the morphological analysis pointed to form design element.

The data statistics stage is also divided into two steps. Firstly, the Kansei word pairs are selected, including the selection of the primary Kansei words and secondary Kansei words that describe the product image. Secondly, the Kansei evaluation of the product image is carried out by the subjects for the experimental samples, and the construction of the product image database is built up through the data statistics. Based on the Kansei Engineering process, the qualitative and quantitative decision models of product image design are established to improve the product image design integrated decision system.

**Figure 4.** The Kansei Engineering experiment process of product image design.

#### **3. Product Image Design Integrated Decision System**

This section may be divided into product image design qualitative decision model, quantitative decision model, integrated decision system and performance evaluation. Firstly, using the Product image design qualitative decision model, the influential design elements for the product image are identified based on Quantification Theory Type I. Secondly, the quantitative decision model is utilized to predict the product total image. GRA–Fuzzy logic sub-models of influential design elements are built up separately. After that, utility optimization model is applied to obtain the multi-objective product image. Thirdly, the product image design integrated decision system is completed to optimize the product image design in the process of product design. Finally, the performance of the system is evaluated by using the root mean square error.

#### *3.1. Product Image Design Qualitative Decision Model*

#### 3.1.1. Quantification Theory Type I

The qualitative decision model of product image design is to identify the main design elements affecting product image by quantitative analysis method. The experimental subjects in Section 2 are divided into four groups. The third group is composed of the professional product designers with rich experience in product design. Through the investigation and statistics of professional designers, the four design elements (product form, color, material and texture) are selected for further research of qualitative decision model.

Quantification Theory Type I is the method that studies the relationship between a set of qualitative variables *x* (independent variables) and a set of quantitative variables *y* (dependent variables) [10]. The mathematical model between variables is built up by multiple regression analysis. According to the mathematical model of Quantification Theory Type I, the procedure for establishing the Kansei-associated model is as follows:

• The mathematical formula of the Kansei-associated model is defined. The product design element is taken as the independent variable *x* (qualitative variable) and the Kansei evaluation value is taken as the dependent variable *y* (quantitative variable). The purpose of Quantification Theory Type I is to solve the regression coefficient of each design element [11]. Therefore, the multiple regression equation of the product design element and Kansei evaluation value can be defined as:

$$\mathfrak{H}\_s^k = \sum\_{i=1}^D \sum\_{j=1}^{M\_i} \beta\_{ij} \mathfrak{x}\_{ijs} \tag{1}$$

where *k* represents the *k*th Kansei image, *k* = 1, 2, ... , *m*, and *m* represents the number of Kansei images; *s* represents the *s*th experimental sample, *s* = 1, 2, ... , *n*, and *n* represents the number of experimental samples; *i* represents the *i*th design element, *i* = 1, 2, ... , *D*, and *D* represents the number of design elements; *j* represents the *j*th category of the *i*th design element, *j* = 1, 2, ... , *Mi*, and *Mi* represents the number of category for the *i*th design element.

*y*ˆ*k <sup>s</sup>* is the evaluation value of the *k*th Kansei image for the *s*th experimental sample. *xijs* represents the independent variable of the *j*th category of the *i*th design element for the *s*th experimental sample, so *xijs* can be defined as follows:

$$\mathbf{x}\_{ijk} = \begin{cases} 1 & \text{the } j \text{th category of the ith design element } f \text{or the } sh \text{ experimental sample} \\ 0 & \text{otherwise} \end{cases} \tag{2}$$

And *xijs* meets the condition of *Mi j*=1 *xijs* = 1, *f or* ∀*i and s*.

β*ij* represents the partial regression coefficient of the *j*th category of the *i*th design element. • Using the least-squares method, error equation can be defined as below:

$$Q = \sum\_{s=1}^{n} \left( y\_s^k - \mathfrak{Y}\_s^k \right)^2 \tag{3}$$

where *Q* is the sum of square error between the predicted value and the actual value. In order to improve the prediction accuracy and minimize the error, the following equation can be defined for the calculation of partial differential:

$$\frac{\partial \mathbb{Q}}{\partial \beta\_{ij}} = 0, \text{ for } \forall \beta\_{ij} \tag{4}$$

β*ij* can be obtained by the simultaneous equations. Namely, the partial regression coefficient of design element can be acquired. The influence degree of each design element on the Kansei image can be obtained by analyzing the numerical value of β*ij*.


Using Quantification Theory Type I, the partial correlation coefficient of the product design element affecting the Kansei image is obtained, and then the decision coefficient of each design element category is acquired. The partial correlation coefficient represents the contribution of product design element to each Kansei image, and decision coefficient demonstrates the precision of model. In the process of product design, if the designers expect to get the higher evaluation value of specific product image, they should give priority to the product design elements with the higher partial correlation coefficient, properly ignore the design elements with the lower partial correlation coefficient, improve design efficiency, save design cost and optimize product image design process.

#### 3.1.2. Product Image Design Qualitative Decision Model

As shown in Figure 5, the qualitative decision model of product image design is constructed. Firstly, the product design elements are analyzed, including product form, product color, product

material, product texture, product pattern and other design elements. The partial correlation coefficient of product element is obtained through Quantification Theory Type I, and the priority ranking of product design elements affecting Kansei image is obtained.

The design elements for the product Kansei image with the greater impact and less impact are analyzed respectively. Designers should pay more attention to the design elements with greater impact, ignore the design elements with less impact, and build a qualitative decision model of product image design. Then, the design elements with greater influence for the profound study are selected, which provides a theoretical basis for the next product image design quantitative decision model.

**Figure 5.** The qualitative decision model of product image design.

#### *3.2. Product Image Design Quantitative Decision Model*

#### 3.2.1. GRA–Fuzzy Logic Model

Grey Relation Analysis (GRA) is a method to determine the relationship (similarity) between two sets of random sequences in grey system. One is a reference sequence (*x*<sup>0</sup> ∈ *X*) and the other is a comparison sequence *X* = {*x*σ|σ = 0, 1, 2, ... , *n*}, and *X* = {*x*σ|σ = 0, 1, 2, ... , *n*} is a set of grey related elements [12]. At a certain moment, the grey correlation degree of two sets of series can be expressed by grey correlation coefficient *r*(*x*0(*k*), *xi*(*k*)), which is defined as below:

$$\begin{array}{l} \left(\mathbf{x}\_{0}(k), \mathbf{x}\_{i}(k)\right) = \min\_{\mathbf{l}} \min\_{\mathbf{l}} \left| \mathbf{x}\_{0}(k) - \mathbf{x}\_{i}(k) \right| + \left. \mathbf{\zeta} \max\_{\mathbf{l}} \mathbf{max} \mathbf{x}\_{k} \right| \mathbf{x}\_{0}(k) - \mathbf{x}\_{i}(k) \Big| /\\ \left| \mathbf{x}\_{0}(k) - \mathbf{x}\_{i}(k) \right| + \left. \mathbf{\zeta} \max\_{\mathbf{l}} \mathbf{max} \mathbf{x}\_{k} \right| \mathbf{x}\_{0}(k) - \mathbf{x}\_{i}(k) \Big| \end{array} \tag{5}$$

*k* = 1, 2, ... , *n*; *i* = 1, 2, ... , *m*

The grey correlation coefficients of different design sub-elements are different. If*r*(*x*0, *xi*) > *r x*0, *xj* , then design sub-element *xi* is more relevant to product image than design sub-element *xj*. Based on the Kansei image database of each design sub-element affecting the product image, the grey correlation coefficient is used to express the relationship between the product image and design sub-elements. By comparing the numerical value of grey correlation coefficient, the sorting of the influence degree of different design sub-elements on product image is obtained. The most influential design sub-elements are selected as input linguistic variables, and the corresponding product image as output linguistic variable of GRA–Fuzzy logic model.

Based on the GRA–Fuzzy theory, combining with the Kansei Engineering experiment, the product design sub-elements with the greater influence are selected first, and the product image of the experimental samples is evaluated by seven-point Likert scale [13]. Then, the influential design sub-elements are taken as the input fuzzy set, and the evaluation values of product image are taken as the output fuzzy set. The triangular fuzzy functions of input and output linguistic variables are obtained. Through the application of fuzzy logic model in MATLAB, the fuzzy rules are constructed, and the center of maximum (CoM) method is used for defuzzification [14]. Finally, the GRA–Fuzzy logic model of product image design is established to determine the specific Kansei image value of product design.

Combined with the product image database in Section 2, the fuzzy rules of GRA–Fuzzy logic model are constructed. In order to reflect the relationship between the multiple design sub-elements and product images involved in the product design process, the "if ... then ... " rules with multiple fuzzy conditions are used as follows:

$$\begin{aligned} \text{If } \mathcal{X}\_1 \text{ is } A\_1 \text{ and } \mathcal{X}\_2 \text{ is } A\_2 \cdots \text{ and } \mathcal{X}\_n \text{ is } A\_n;\\ \text{Then } \mathcal{Y}\_1 \text{ is } B\_1 \text{ and } \mathcal{Y}\_2 \text{ is } B\_2 \cdots \text{ and } \mathcal{Y}\_n \text{ is } B\_n \end{aligned} \tag{6}$$

where *A*1, *A*2, ... , *An* and *B*1, *B*2, ... , *Bn* are fuzzy linguistic terms, which are, respectively, represented by input linguistic variables *X*1, *X*2, ... , *Xn* and output linguistic variables *Y*1,*Y*2, ... ,*Yn*. The input linguistic variables *X*1, *X*2, ... , *Xn* represent sub-elements of product design, and the output linguistic variables *Y*1,*Y*2, ... ,*Yn* represent product image. Each fuzzy rule of GRA–Fuzzy logic model of product image design is to associate the given combination of product design sub-elements with the corresponding numerical state of product image.

Defuzzification is the process of transforming the membership of output linguistic variables into numerical value. The most commonly used defuzzification technology, the center of maximum (CoM), is determined as follows:

$$\text{yCoM} = \frac{\sum\_{i} \left[ \mu(y\_i \times y\_i) \right]}{\sum\_{i} \mu(y\_i)} \tag{7}$$

where *i* represents the linguistic item of output linguistic variable, *yi* is the maximum value of each linguistic item *i*, and μ(*yi*) is the aggregate output membership function.

Finally, the GRA–Fuzzy logic model of product image design is constructed to assist designers to evaluate the product image design. When the product designers input the combination of product sub-elements, they can output the value of product image and obtain the numerical states of multiple product images.

#### 3.2.2. Utility Optimization Model

In the product design process, it is difficult to judge design scheme by using a single index. For the solution method of multi-objective programming, the utility optimization model is generally used through transforming the multi-objective model into the single-objective model. The model solution basic idea is that each objective function of the programming problems can be calculated in a certain way [15].

$$\max Z = \psi(X) \tag{8}$$

$$\text{s.t.}\,\Phi(X) \le G \tag{9}$$

ψ is the sum function of utility functions related to each objective function.

When using utility function as planning objective, it is necessary to determine a set of *wi* that reflects the weight of each objective function in the overall objective of the original problem, namely:

$$\max \psi = \sum\_{i=1}^{k} w\_i \psi\_i \tag{10}$$

$$\Phi(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_n) \le \mathbf{g}\_i(i = 1, 2, \dots, m) \tag{11}$$

where *wi* should satisfy - *k i wi* = 1

Parametric vector form is shown as below:

$$\mathbf{m} \mathbf{a} \mathbf{x} \psi = \mathbf{w}^T \psi \tag{12}$$

$$\text{s.t.}\,\Phi(X) \le G \tag{13}$$

#### 3.2.3. Product Image Design Quantitative Decision Model

The quantitative decision model of product image design is constructed as shown in Figure 6. Based on the Kansei Engineering experiment, the qualitative decision model of product image design in Section 3.1.2 is used to get the influential design elements for product image, and the GRA–Fuzzy logic sub-models of the influential design elements are constructed respectively.

**Figure 6.** The quantitative decision model of product image design.

Firstly, using the Kansei Engineering experiment process, the design elements with great influence are analyzed, and the sub-elements of the influential design elements are obtained. Meanwhile, the representative Kansei words of the product image expected by users are determined. The subjects evaluate the product image of the experimental samples by Likert scale to set up the product image database. Then, through the grey relational analysis (GRA), the grey correlation coefficients of the design sub-elements that affecting the product image are obtained, and the numerical values of the grey correlation coefficients are sorted and compared in order to obtain the design sub-elements that have a great influence on the product image.

The design sub-elements with great influence obtained above are taken as input linguistic variables, and the product image expected by users is taken as an output linguistic variable. The fuzzy functions of input linguistic variables and output linguistic variables are constructed by the function of triangular membership, and the GRA–Fuzzy logic model is constructed by MATLAB. The fuzzy control data

table and the corresponding fuzzy rules between design sub-elements and product image are obtained. Then, the fuzzy numbers are defuzzified by using the CoM method. Finally, the GRA–Fuzzy logic sub-models of the influential design elements are established.

Users usually expect to obtain multiple product images in the practical design and production process. The utility optimization model is used to transform the multi-objective product image into the single-objective product image, optimize the GRA–Fuzzy logic model of product design elements. The weights of product image are obtained according to the demand of design objectives, which reflects the weight of each function of product image in total product image objective, and the product total image value is obtained by linear weighting method. The quantitative decision model of product image design is optimized, namely:

$$\log(\mathbf{x}) = \sum\_{j=1}^{q} w\_j \mathbf{g}\_j(\mathbf{x}) \tag{14}$$

where *gj*(*x*) is the GRA–Fuzzy logic model of product image *j*. *wj* is the weight of the *j*th product image function *gj*(*x*) obtained according to the requirements of objective programming in the product total image objective.

The GRA–Fuzzy logic sub-models of product design elements can help designers evaluate the design elements of product image. Based on the construction of GRA–Fuzzy logic model and utility optimization model, the product image design model is further optimized, and the quantitative decision model of product image design is constructed.

#### *3.3. Product Image Design Integrated Decision System*

As shown in Figure 7, product image design integrated decision system is composed of qualitative decision model and quantitative decision model. For the varied kinds of industrial product design, firstly, the qualitative decision model of product image design is constructed to obtain the product design elements that have a greater influence on the user's Kansei image. After that, the quantitative decision model of product image design is implemented to acquire the relationship between the influential design sub-elements and the numerical value of Kansei image. The collaborative work of the above two models constitutes the integrated decision system of product image design. Therefore, it is more efficient and accurate that the combination of product design elements matches with the output of user's Kansei image.

**Figure 7.** The product image design integrated decision system.

#### *3.4. Performance Evaluation*

In order to evaluate the performance of the product image design integrated decision system, the root mean square error (*RMSE*) is used as follows [16]:

$$RMSE = \sqrt{\frac{\sum\_{i=1}^{n} \left(\mathbf{x}\_{i} - \mathbf{x}\_{0}\right)^{2}}{n}} \tag{15}$$

where *xi* is the *i* th output value of product image design integrated decision system and *x*<sup>0</sup> is the expected value evaluated by the experimental subjects. We use test samples to evaluate the prediction performance of product image design integrated decision system, and the subjects of the second group evaluate the product total image of test samples to get the numerical value. Then, the design sub-elements of each sample are input into the product image design integrated decision system, and the output value of the product total image is obtained, which compared with the evaluation value of the subjects. The *RMSE* of product image design integrated decision system is obtained. Then, the performance evaluation of product image design integrated decision system is completed.

The case study shows that the optimized product image design integrated decision system has a good performance for predicting and evaluating the product image. The proposed system has a better prediction ability of Kansei image for the new product with given design elements.

#### **4. Case Study of Train Seat Image Design**

Based on the research of the optimization process of the product image design integrated decision system in Section 3, this section aims at optimization of the high-speed train seat image design decision system. The qualitative and quantitative decision models of the high-speed train seat image design are constructed to further promote the process of the train seat image design integrated decision system.

#### *4.1. Kansei Engineering Experiment*

#### 4.1.1. Design Investigation Stage

Firstly, the experimental subjects and samples are selected. Afterwards, the final experimental samples are determined by the selection of train seat form, color, material and texture. Finally, the representative Kansei image words for describing the train seat image are chosen.

• Experimental subject selection

The experimental study involves 60 subjects, divided into four groups. The first group includes general users (10 males and 10 females). In the second group, there are expert users (10 males and 10 females). The third group is composed of 10 professional product designers (6 males and 4 females), and the designers of this group have at least 6 years of product design experience on average. The fourth group includes 10 experts (5 males and 5 females) in the field of high-speed train.

• Experimental sample selection

The subjects of the first group classify the selected product pictures based on their similarity degree of appearance design of train seat, using the Kawakita Jiro (KJ) method. Therefore, the 20 subjects of the first group classify 198 seat pictures of different manufacturers and models that entered the market from 2009 to 2019. Based on the classification results of the subjects, 30 representative experimental samples of train seats are selected by multi-dimensional scale analysis and cluster analysis.

• Design element analysis

The train seat form samples are extracted from 30 representative seat samples. According to the professional knowledge and design experience, 10 designers of the third group obtain the four design elements (including the seat form, color, material and texture) for future analysis. After that, they are required to extract seat form samples through similarity classification, and then 10 subjects are organized into focus groups to summarize and integrate their research results. After that, the 6 form samples are obtained, which are divided into three groups on average, the first group is the form experimental samples. The form experimental sample of train seat No. 2 of the first group is shown in Figure 8.

**Figure 8.** The form experimental sample of high-speed train No. 2 of the first group.

This study regards a small amount of color as a huge work of color speculation and establishes the relationship between product form, color, material, texture and product image [17]. The 10 subjects are organized into a focus group to conduct the inductive integration and analysis of color selection. Finally, the 6 color samples are selected as the color experimental samples of train seat. The form samples are filled with the 6 color experimental samples respectively. The form sample No. 2 is filled with the 6 color experimental samples in Figure 9 as below.

**Figure 9.** The 6 color experimental samples of the train seat form sample No. 2.

The subjects of third group analyzed the material design and texture design of the train seat, and the representative material samples and texture samples of train seats are collected and selected by the third subjects, and the material samples and texture samples are applied to the color samples. Through the research and investigation, the representative material samples are divided into cotton fabric, wool fabric, chemical fiber and leather, as shown in Table 1.


The representative texture samples are divided into coarse texture and fine texture. The above 4 material samples and 2 texture samples are applied to 12 train seat samples in the first group, a total of 96 train seat experimental samples, as the final experimental samples of this study.

#### • Kansei image word selection

Based on the preceding steps, three groups of representative Kansei words are determined to describe the Kansei image of train seat: traditional–modern, simple–complex, unique–ordinary. We use "traditional–modern" (T–M) as the primary Kansei word of the train seat, which is the target image for constructing the qualitative decision model and quantitative decision model of the seat image design. The above experimental process becomes the application case for developing the method of Kansei Engineering.

#### 4.1.2. Data Statistics Stage

The expert users evaluate the product image matched with the seat samples. As shown in Figure 10, semantic difference method is used to make the seven-point Likert scale, and 20 subjects are required to carry out the numerical value evaluation of T–M product image for the experimental samples (1 and 7 represent the extremely traditional and extremely modern image respectively). For example, seat No. 2 has "modern" image, namely that the maximum value of T–M image is 7, the minimum value is 3, and the average value is 5.33. The data statistics stage provides the experimental database for the qualitative decision model and quantitative decision model of seat image design.



#### *4.2. Train Seat Image Design Qualitative Decision Model*

Through the analysis of Quantification Theory Type I, the importance and influence degree of train seat design elements are obtained. The independent variables (such as seat form, color, material and texture) of factors quantification are the nominal scales, so the numerical values of 1, 2, 3 and 4 are used to represent the categories of design elements, for example, "cotton fabric", "wool fabric", "chemical fiber" and "leather" of material samples are, respectively, represented by 1, 2, 3 and 4. The result is shown in Table 2.


**Table 2.** Analysis of result from Quantification Theory Type I.

The partial correlation coefficients of design elements in Table 2 show that the correlation between four independent variables (seat form, color, material and texture) and Kansei image words (T–M). For example, for the T–M image, the partial correlation coefficients of design elements are 0.441, 0.791, 0.290 and 0.003, respectively, X2 of which is the highest, indicating that the "color" factor is the design factor that has the greatest influence on T–M image of the independent variables. The second factor is seat form (X1), which shows that when the design expectation of the train seat is to obtain the T–M image, the designer should pay more attention to color design element and form design element. On the contrary, designers can properly ignore the design elements with less relevance, such as seat material (X3) and texture (X4), because these two design elements have a relatively small impact on T–M image.

Based on the above experimental process of Kansei Engineering, the qualitative decision model of train seat image design is finally obtained by using Quantification Theory Type I method. The qualitative decision model of train seat image design is shown in Figure 11. The partial correlation coefficients of seat design elements are obtained through Quantification Theory Type I. The sequence of design elements affecting the Kansei T–M image of train seat is color, form, material and texture. Therefore, the study selects the design elements that have a great influence on the Kansei image: seat form and seat color, which provides a theoretical basis for constructing the quantitative decision model of train seat image design in the next step.

**Figure 11.** The qualitative decision model of train seat image design.

#### *4.3. Train Seat Image Design Quantitative Decision Model*

Based on the conclusion of train seat image design qualitative decision model, the design elements that have great influence on seat image are seat form and color design elements. Using the experimental process of Kansei Engineering, the GRA–Fuzzy logic sub-models of train seat form design and color design are respectively constructed. Based on the GRA–Fuzzy theory, this paper analyzes the form sub-elements and the color sub-elements. Aimed at the representative Kansei words of the seat image, the subjects evaluate the Kansei image of seat samples through the seven-point Likert scale to build up the train seat image database.

As the result of design element analysis, the 9 form sub-elements and their corresponding element types are extracted from the 30 representative train seat samples through morphological analysis. In this section, we select one representative train seat randomly from the 30 representative train seat samples and fill the selected train seat with 72 representative colors as the experimental train seat color samples, for a total of 72 color samples. In Table 3, Column 1 and 7 display the numbers of color samples, and the remaining columns display the corresponding color sub-element ("Hue (*h*◦)" element (*xc*1), "Chroma (*C*\*)" element (*xc*2), "Lightness (*L*\*)" element (*xc*3), "parameter *a*\*" element (*xc*4) and "parameter *b*\*" element (*xc*5)) values of these color samples.


**Table 3.** Color element values of the representative train seat color sample C22–C33.

The GRA method is used to obtain the grey correlation degree of form design sub-element as shown in Table 4 from Section 3.2.1, and the sub-elements of form design that affect the Kansei image of train seat are, respectively, "Seat back waistline" (*x*5), "Seat back body shape" (*x*4), "Headrest shape" (*x*2), "Length ratio between seat-back to seat" (*x*8), "Relationship between seat-back and headrest" (*x*3) and "Seat-back slope angle" (*x*9). The above form design sub-elements that have greater influence on Kansei image are taken as the input linguistic variables and the Kansei image of train seat as output linguistic variables. The GRA–Fuzzy logic sub-model between the form design sub-elements and Kansei image is constructed by the triangular membership function.


**Table 4.** The numerical values of *r*(*x*0, *xi*).

Afterwards, the corresponding 54 fuzzy rules are built up, and the GRA–Fuzzy logic sub-model of train seat form design is constructed. The fuzzy rules No. 22–No. 30 have been added as followed in Table 5.

**Table 5.** The fuzzy rules No.22–No.30 of GRA–Fuzzy logic sub-model of form design sub-elements.


According to Section 3.2.1, GRA method is used to obtain the grey correlation degree of color design sub-element. It shows that the color design sub-elements that affect the Kansei image of train seat are "Lightness (*L*\*)" (*xc*3), "Chroma (*C*\*)" (*xc*2), "parameter *a*\*" (*xc*4) and "parameter *b*\*" (*xc*5). Taking the above color design sub-elements that have great influence on seat Kansei image as input linguistic variables and seat Kansei image as output linguistic variables, the GRA–Fuzzy logic sub-model of color design sub-elements matched with Kansei image is constructed by triangular membership function. Next, the 144 corresponding fuzzy rules are established, and the GRA–Fuzzy logic sub-model of train seat color design is constructed.

Based on the cooperative work of GRA–Fuzzy logic sub-models of the aforementioned train seat form design and color design, in terms of the application of utility optimization model, the utility optimization model is constructed for the Kansei image output value of train seat for perceptual cognition of users. Three groups of representative Kansei words are selected, including traditional–modern (T–M), simple–complex (S–C) and unique–ordinary (U–O).

In accordance with the method of multi-image objective programming, three corresponding GRA–Fuzzy logic sub-models of T–M, S–C and U–O images are established, respectively. When using utility function as the planning objective, three corresponding weights *w*1, *w*<sup>2</sup> and *w*<sup>3</sup> are used to reflect the weight of each GRA–Fuzzy logic sub-model of each Kansei image in the overall objective, and the utility optimization model of the product total image is obtained. When the primary image of seat design expected by users is T–M image, followed by S–C and U–O image, according to the planning of users' perceptual needs, the weight of GRA–Fuzzy logic sub-model for T-M image is 0.4, and the weight of GRA–Fuzzy logic sub-model reflecting S–C image and U–O image is 0.3, respectively. Finally, the seat design total image of multi-image planning objective of train seat is obtained.

Based on the experiment process of Kansei Engineering, combined with the utility optimization model, the seat image design quantitative decision model of high-speed train is finally constructed, and the specific process is shown in Figure 12.

**Figure 12.** The quantitative decision model of train seat image design.

#### *4.4. Train Seat Image Design Integrated Decision System*

As shown in Figure 13, the train seat image design integrated decision system consists of qualitative decision model and quantitative decision model of seat image design. Aiming at the seat image design of high-speed train, the qualitative decision model of the seat image design is constructed firstly, and then the seat form design elements and color design elements that affect the Kansei image are obtained. After that, the quantitative decision model of the seat image design is constructed to obtain the predicted numerical value of the total image of the corresponding seat. The collaborative work of the above two models constitutes the integrated decision system of train seat image design, which makes the image design process of train seat more intelligent, accurate and efficient.

**Figure 13.** The train seat image design integrated decision system.

#### *4.5. Performance Evaluation*

In order to verify the prediction ability of train seat image design integrated decision system, five test samples shown in Table 6 are used, including form design sub-elements and color design sub-elements that have greater influence on product image. The performance of train seat image design integrated decision system is evaluated by using the root mean square error.


**Table 6.** Input and output values of the seat test samples.

We select T–M image as the primary image of seat design, then S–C and U–O image as the secondary image, and apply the seat image design integrated decision system to predict the product total image. The second row of Table 7 shows the average value of the product total image evaluated by the 20 subjects of second group on five test samples, and meanwhile, the third row indicates the prediction value of the integrated decision system of seat image design. The results show that the seat image design integrated decision system has lower *RMSE* value. Therefore, the integrated decision system of seat image design provides an effective mechanism for train seat design to match with the total image of product design expected by users.

**Table 7.** *RMSE* result of the integrated decision system of train seat total image design.


The good performance of the seat image design integrated decision system shows that the system can help the designer to evaluate the combination of train seat design elements, aiming at satisfying the expected Kansei image. Therefore, this paper provides an effective design evaluation system for the improvement of seat design reflecting the ideal product image through the product image design integrated decision system. For example, the product designer can input the values of multiple design sub-elements by using the product image design integrated decision system, and then the predicted value of multi-objective product image is obtained. If the predicted value does not meet the standard of designer, the corresponding new product image predicted value can be obtained by modifying the combination of design sub-elements until the designer obtains the satisfactory product image value.

#### *4.6. Computer-Aided Seat Design and Image Evaluation*

Combined with the computer-aided seat form design and color design, the seat form and color design elements are input into the train seat image design integrated decision system, and then the product image value is obtained. Through the further application of the computer-aided product form modeling and color rendering, train seat image design integrated decision system is verified for the predictable performance of product image.

Based on the computer-aided seat test sample form design and color design, four color design samples are applied to three form design samples, and a total of 12 product design samples are obtained. Figure 14 shows the seat rendering of the form sample 3 filled with the color sample 1.

**Figure 14.** The train seat rendering of the form sample 3 filled with the color sample 1.

At first, the experimental subjects are required to evaluate the T–M image of the above 3-D seat models. After that, the T–M image outputs of the seat models are carried out by using the train seat image design integrated decision system. The T-M image results of models 5 to 10 are shown in the third row of Table 8. Finally, the *RMSE* value of the train seat image design integrated decision system is calculated with the result of less than 0.5, which proves that integrated decision system of train seat image design performs better in predicting the T–M image.

**Table 8.** *RMSE* result of the integrated decision system of train seat T–M image design.


In conclusion, the integrated decision system of seat image design can effectively predict the image value of 3-D seat model of computer-aided design and rendering. Therefore, utilizing computer software to build 3-D model of industrial product, combined with computer-aided 3-D printing rapid prototyping, we can assist the use of product image design integrated decision system, which makes the image evaluation of product design more intuitive and efficient. Figure 15 illustrates the 3-D printed drawing of the form sample 3. In this study, 3-D model printing is carried out according to the ratio of 5:1 reduction of seat 3-D model to actual size in order to assist the evaluation process of integrated decision system of seat image design.

**Figure 15.** The 3-D printed drawing of the seat form sample 3.

#### **5. Research Limitations and Discussions**

This paper takes the seat image design of high-speed train as the research case, and this method is also applicable to other industrial product design, interaction design, service design and user research fields. Based on the experimental framework of Kansei Engineering, the Kansei image perceived by users is transformed into product design elements, interaction design elements or service design elements, etc. Moreover, by establishing quantitative and qualitative decision models of product image design, the integrated decision system of product image design is obtained. The research system can evaluate the perceptual experience of users, make the designers work more efficiently and accurately, and promote and optimize the process of product design, interaction design, service design and user research by predicting the Kansei image of products.

The next work is to further explore the product image design integrated decision system based on Kansei Engineering, combined with product environment elements. The integrated decision system will be expanded and upgraded in the future. With the change of society and the passing of time, the perceptual cognition of users will change constantly. The evaluation database of experimental subjects and the Kansei image database of products in the integrated decision system of product image design need to be adjusted and updated timely to ensure that the integrated decision system of product image design can go with the tide of development and the user environment for any new product market.

In the intelligent application environment, more high and new technologies are developing rapidly. Therefore, the development trend of the product image design integrated decision system will be more intelligent, accurate, digital and systematic. Meanwhile, in the establishment process of experimental database, the design concept, 3-D model and rendering of product design will be further optimized and promoted [18]. Therefore, the intelligent research of product image design integrated decision system has broad prospects and far-reaching research significance.

#### **6. Conclusions**

In this paper, we have proposed a methodology of product image design integrated decision system based on Kansei Engineering theory. Firstly, Quantification Theory Type I is used to effectively identify the influential design elements that affect the product image, and the product image design qualitative decision model is constructed. Secondly, according to GRA–Fuzzy theory, the GRA–Fuzzy logic sub-models of influential design elements for the product image are constructed. Then, combining with the utility optimization model, the product image design quantitative decision model is built up. Finally, the product image design integrated decision system is constituted to analyze, decide and evaluate the product image design. The result shows that the system is effective in predicting the product images and optimizes the product image design process.

To illustrate the product image design integrated decision system, we have performed an experimental study on train seat design. The result demonstrates that the system has a better performance for improving the product image design. The designers can effectively plan the product design for the specific product image to meet the user perceptual expectation. Although the train seat design is used as a case in this paper, the methodology can be applied to other user-oriented product designs with multiple design elements. Along with the development of modern science and technology innovation and development, the product image design integrated decision system involved with multiple product design elements and product images can be further improved and developed in the future.

**Author Contributions:** Conceptualization, L.X. and X.Y.; methodology, L.X.; software, L.X.; validation, L.X. and Y.Z.; formal analysis, L.X.; investigation, Y.Z.; resources, X.Y.; writing—original draft preparation, L.X.; writing—review and editing, L.X., X.Y. and Y.Z.; visualization, L.X. and Y.Z.; supervision, X.Y.; project administration, X.Y.; funding acquisition, L.X. and X.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Major Program of Beijing Municipal Science and Technology Commission under grant number D16111000110000.

**Acknowledgments:** We would like to thank the 60 subjects, 10 product designers, 10 design experts for their participation and assistance in the experimental study. We also appreciate the editor and anonymous reviewers for their valuable comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **A Coordination Space Model for Assemblability Analysis and Optimization during Measurement-Assisted Large-Scale Assembly**

#### **Zhizhuo Cui and Fuzhou Du \***

School of Mechanical Engineering & Automation, Beihang University, Beijing 100191, China; by1407136@buaa.edu.cn

**\*** Correspondence: du\_fuzhou@buaa.edu.cn; Tel.: +86-1352-182-3595

Received: 17 April 2020; Accepted: 8 May 2020; Published: 11 May 2020

#### **Featured Application: For assembly incoordination caused by excessive assembly deviations, the proposed method can predict the assemblability and solve the assembly features that need accuracy compensation, to improve the assembly e**ffi**ciency.**

**Abstract:** The assembly process is sometimes blocked due to excessive dimension deviations during large-scale assembly. It is inefficient to improve the assembly quality by trial assembly, inspection, and accuracy compensation in the case of excessive deviations. Therefore, assemblability prediction by analyzing the measurement data, assembly accuracy requirements, and the pose of parts is an effective way to discover the assembly deviations in advance for measurement-assisted assembly. In this paper, a coordination space model is constructed based on a small displacement torsor and assembly accuracy requirements. An assemblability analysis method is proposed to check whether the assembly can be executed directly. Aiming at the incoordination problem, an assemblability optimization method based on the union coordination space is proposed. Finally, taking the space manipulator assembly as an example, the result shows that the proposed method can improve assemblability with a better assembly quality and less workload compared to the least-squares method.

**Keywords:** measurement-assisted assembly; coordination space; assemblability; small displacement torsor

#### **1. Introduction**

Large-scale mechanical products like ships, automobiles, aircrafts, etc. are complex in structure, large in size, and accurate in assembly quality. The assembly workload of the manufacturing process is heavy [1]. These products often need accuracy compensation in the assembly process because of the excessive assembly deviations, which lead to inefficiency. The assembly deviations might be caused by the eventual poor machining quality of parts, or excessive tolerances set by designers. Thus, the trial assembly is often used to detect the assembly deviations in advance, and the parts are then separated to make an accuracy compensation on the bad dimensions. The assembly process takes a long time by the following steps: Trial assembly, measurement of deviations, separation of parts, and re-trial assembly. Therefore, an assemblability analysis and optimization method based on the measurement data is necessary to predict the assembly deviation and make the accuracy compensation in advance.

With the development of measurement-assisted assembly (MAA) [2], measurement technology has become a bridge between the real world and the digital world. Marguet et al. [3] introduced a MAA application in an airbus assembly line. The least-squares method was used to calculate the optimal pose. Chen et al. [4] proposed a weighted SVD algorithm to obtain the optimal pose of components, which improved the accuracy of pose evaluation. Li et al. [5] proposed a coaxial alignment method

using distributed monocular vision. The iterative reweighted particle swarm optimization method was constructed to improve the measurement ability of complicated wearing holes. Wang et al. [6] calculated the assembly clearance of a wing-fuselage assembly based on the optimal pose. The above methods mainly consider the measurement and calculation of the assembly pose, and then realize alignment through pose adjustment tooling. The assembly will be difficult if the quality of the parts is poor.

Assemblability prediction is the first step to judge whether the assembly is qualified in the measurement-assisted assembly. Sukhan et al. [7] evaluated the assemblability based on tolerance propagation. Sanderson et al. [8] assessed the assemblability by the maximum likelihood problem, which was solved by the Kalman filter algorithm. The traditional assemblability evaluation methods are mainly used to find the assembly problem in the design phase, but not in the assembly phase. Cui and Du [9] proposed the concept of pose feasible space to assess the assembly coordination. Yuan et al. [10] proposed an assembly quality assessment method based on weighted geometric constraints to calculate the optimal pose. Wu et al. [11] proposed a constraint coordination index to assess the assembly quality. Ma et al. [12] developed the assembly precision pre-analysis technique in the simulation of virtual assembly. Du et al. [13] proposed a pose decoupling model of the axis tolerance feature to decouple the analysis of any pose within the tolerance domain.

The accuracy compensation methods are used to improve assemblability. The digital compensation method has become a research highlight to improve the assemblability. Davis et al. [14] put forward the method of measuring the assembly clearance and realizing the digital manufacturing of the accuracy compensation gasket. Fabian et al. [15] introduced a shimming method by 3D printing technology, and the assembly clearance was measured by optical measurement. Wang et al. [16] provided a shimming method based on scanned data for a wing box assembly involving non-uniform gaps. In addition, finite element analysis was taken to improve the shimming scheme. Those methods, however, need to be assembled first, followed by measurement of the deviations to be compensated, resulting in a lower efficiency.

Some scholars proposed predictive shimming and predictive fettling methods to improve the assembly efficiency and quality [17]. Cui et al. [18] proposed the oriented points group to calculate the deviation of multiple shaft-and-holes, and the gap was shimmed. Yang et al. [19] analyzed the deviation from the measured point cloud to the model to improve skin finishing. Yu et al. [20] employed a virtual assembly and repair analysis method based on both the geometric design model and object scanning model. Manohar et al. [21] proposed an alternative strategy for predictive shimming, based on machine learning and sparse sensing to first learn gap distributions from historical data. Lei et al. [22] presented an automated and in situ alignment approach with the assistance of computer numerical controlled (CNC) positioners and laser trackers to reduce the finish machining workload. The above studies are aimed at specific cases.

The accuracy compensation method is usually applied after assembly. Then, the assembly sometimes needs be separated, which leads to low efficiency. In this paper, an assemblability analysis and optimization method based on the coordination space model is constructed during measurement-assisted large-scale assembly. In Section 2, the coordination space model based on the small displacement torsor is constructed. In Section 3, the assemblability analysis based on the coordination space model is proposed. In addition, the uncoordinated case is further analyzed. In Section 4, the assemblability optimization method based on the union coordination space is proposed for the uncoordinated case. In Section 5, the space manipulator assembly is taken as an example to verify the proposed method. The result shows that the proposed method can optimize the assemblability with less workload and better assembly quality compared to the least-squares method.

#### **2. Coordination Space Model Based on Small Displacement Torsor**

Assemblability refers to the ability of parts to satisfy the assembly accuracy requirements in terms of dimensions, which can be expressed by coordination accuracy. Traditionally, coordination accuracy [23] is the difference in the manufacturing dimensions. Figure 1 shows the coordination accuracy of a keyway assembly.

**Figure 1.** Coordination accuracy of a keyway assembly.

The coordination accuracy is

$$
\nabla\_{AB} = \mathbf{1}\_1 + \mathbf{2} = L\_A - L\_{B'} \tag{1}
$$

It can be seen that the coordination accuracy is the amount of the allowance on a certain dimension. In this way, the assembly coordination of a single dimension is well presented by coordination accuracy such as angle, length, etc. However, it is not suitable for complicated assembly. Therefore, the concept should be extended to pose allowance space and the space can be predicted by digital measurement data during large-scale assembly. This space is named the assembly coordination space, which is the ability of pose variation under the condition of assembly accuracy requirements, as Figure 2 shows.

**Figure 2.** The pose variations under the condition of assembly accuracy requirements.

The parts of the assembly are divided into the reference part and the align part. The reference part is the fixed part during assembly and the align part will move to the target pose by the pose adjustment tooling. Assume that the primary measurement data of the two parts are

$$\begin{cases} \begin{aligned} \mathbf{P}^{\mathbb{R}} &= \begin{bmatrix} p\_1^{\mathbb{R}} \ p\_2^{\mathbb{R}} \dots \ p\_n^{\mathbb{R}} \end{bmatrix} \\ \mathbf{P}^{A} &= \begin{bmatrix} p\_1^{A} \ p\_2^{A} \dots \ p\_n^{A} \end{bmatrix} \end{aligned} \tag{2}$$

where *P<sup>R</sup>* and *P<sup>A</sup>* are the point sets of the reference part and align part, separately, where *p<sup>R</sup>* <sup>1</sup> , etc. and *pA* <sup>1</sup> , etc. are the points of the sets *<sup>P</sup><sup>R</sup>* and *<sup>P</sup>A*, respectively. The two parts are separated first. According to the least-squares method, the optimal assembly pose can be calculated by

$$\begin{cases} \begin{aligned} \mathbf{P}^{\mathbb{R}} & \approx \mathbf{R} \mathbf{P}^{A} + T\\ \boldsymbol{\varepsilon} &= \min \left\{ \sum\_{i=1}^{n} ||\mathbf{R} \mathbf{p}\_{1}^{A} + T - \mathbf{P}^{\mathbb{R}}||^{2} \right\} \end{aligned} \tag{3}$$

where *R* is the rotation matrix, *T* is the movement matrix, and *e* is the minimum residual sum of squares. The singular value decomposition method [24] is taken to calculate the parameter *R* and *T*. Then, the optimal pose of the align part based on the least-squares method is

$$
\omega\_0 = \begin{bmatrix} \mathbf{R} & T \\ 0 & 1 \end{bmatrix} \tag{4}
$$

The assembly deviation can be predicted by the pose ω<sup>0</sup> of the align part. The key assembly characteristics (KAC) [25] are the important geometric structures that have key influences on assembly quality. They are described by measurement data and some dimensions that are not necessary to be measured.

$$\mathbf{K} = \langle P, \mathbf{G} \rangle\_{\prime} \tag{5}$$

where *K* is the parameters of a KAC, *P* is the measurement data, *G* is the dimensions that are not necessary to be measured. The KACs have an irregular distribution in space during large-scale assembly. As shown in Figure 3, the wing-fuselage assembly is completed by 4 pairs of joints. There are four assembly accuracy requirements on each pair of joints: Two on coaxialities and two on clearances.

**Figure 3.** Wing-fuselage assembly.

The KACs are restrained by assembly accuracy requirements. The assembly accuracy is described as

$$\mathbf{^i T^i\_j} = f^i\_j(\mathbf{K^R\_i}, \mathbf{K^A\_i})\_\prime \tag{6}$$

where *K<sup>R</sup> <sup>i</sup>* is the parameters of the *<sup>i</sup>*th KAC on the reference part (fuselage), *<sup>K</sup><sup>A</sup> <sup>i</sup>* is the parameters of the *i*th KAC on the align part (wing), *T<sup>i</sup> <sup>j</sup>* is the *<sup>j</sup>*th assembly accuracy of the *<sup>i</sup>*th KAC, and *<sup>f</sup><sup>i</sup> <sup>j</sup>* is the mapping from parameters to the assembly accuracy. The assembly accuracy should meet the requirements of assembly accuracy, which is formulated in Equation (7)

$$T\_j^i \in \left[ T\_j^{i-\min}, T\_j^{i-\max} \right],\tag{7}$$

where *Ti*−*min <sup>j</sup>* and *<sup>T</sup>i*−*max <sup>j</sup>* are the ranges of *<sup>T</sup><sup>i</sup> j* . Substitute Equation (6) into Equation (7):

$$\{f\_j^i(\mathbf{K}\_i^R, \mathbf{K}\_i^A) \in \left[T\_j^{i-\min}, T\_j^{i-\max}\right]\_\prime\}\tag{8}$$

For the m assembly accuracy requirements on n KACs of the assembly, the constraint equations can be expressed as

$$\forall f \big( \mathbf{K}\_i^R, \mathbf{K}\_i^A \big) \in \left[ T\_j^{i-\min}, T\_j^{i-\max} \right], i \in [1, n], j \in [1, j\_i], \tag{9}$$

where *ji* is the assembly accuracy requirement number of the *i*th KAC, *m* = *n <sup>i</sup>*=<sup>1</sup> *ji*. When all KACs satisfy their assembly accuracy requirements, the pose is a valid pose to be aligned.

As shown in Figure 4, the valid pose may not be the only one that satisfies all assembly accuracy requirements. Therefore, the adjacent poses of the primary pose shown in Figure 4a can be analyzed. A small displacement torsor (SDT) [26] represents a tiny rigid body's pose variation. It is described as

$$
\omega\_{\Delta} = (x\_{\Delta \prime} y\_{\Delta \prime} z\_{\Delta \prime} \alpha\_{\Delta \prime} \beta\_{\Delta \prime} \gamma\_{\Delta}),
\tag{10}
$$

**Figure 4.** Valid poses of wing-fuselage assembly.

The homogeneous transformation matrix of an SDT is

ω*<sup>H</sup>* <sup>Δ</sup> = ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ *<sup>C</sup>*γΔ*C*β<sup>Δ</sup> <sup>−</sup>*S*γΔ*C*β<sup>Δ</sup> *<sup>S</sup>*β<sup>Δ</sup> *<sup>x</sup>*<sup>Δ</sup> *<sup>S</sup>*γΔ*C*α<sup>Δ</sup> <sup>+</sup> *<sup>C</sup>*γΔ*S*βΔ*S*α<sup>Δ</sup> *<sup>C</sup>*γΔ*C*α<sup>Δ</sup> <sup>−</sup> *<sup>S</sup>*γΔ*S*βΔ*S*α<sup>Δ</sup> <sup>−</sup>*C*βΔ*S*α<sup>Δ</sup> *<sup>y</sup>*<sup>Δ</sup> *<sup>S</sup>*γΔ*S*α<sup>Δ</sup> <sup>−</sup> *<sup>C</sup>*γΔ*S*βΔ*C*α<sup>Δ</sup> *<sup>C</sup>*γΔ*S*α<sup>Δ</sup> <sup>+</sup> *<sup>S</sup>*γΔ*S*βΔ*C*α<sup>Δ</sup> *<sup>C</sup>*βΔ*C*α<sup>Δ</sup> *<sup>z</sup>*<sup>Δ</sup> 0 0 01 ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ ≈ ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ 1 −γ<sup>Δ</sup> β<sup>Δ</sup> *x*<sup>Δ</sup> γ<sup>Δ</sup> 1 −α<sup>Δ</sup> *y*<sup>Δ</sup> −β<sup>Δ</sup> α<sup>Δ</sup> 1 *z*<sup>Δ</sup> 0 0 01 ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ , (11)

where *S* is sin, *C* is cos, lim αΔ→0 *C*α<sup>Δ</sup> = 1, lim αΔ→0 *<sup>S</sup>*α<sup>Δ</sup> <sup>=</sup> <sup>α</sup>Δ, and lim <sup>α</sup>Δ,βΔ→<sup>0</sup> *S*αΔ*S*β<sup>Δ</sup> = 0. The change in point *p* = (*x*, *y*, *z*) after a slight change in the rigid body's pose is

$$p'\_{cx} = [\mathbf{x} \ \mathbf{y} \ z \ \mathbf{1}] \begin{bmatrix} 1 & -\boldsymbol{\gamma}\_{\Lambda} & \boldsymbol{\beta}\_{\Lambda} & \mathbf{x}\_{\Lambda} \\ \boldsymbol{\gamma}\_{\Lambda} & 1 & -\boldsymbol{\alpha}\_{\Lambda} & \mathbf{y}\_{\Lambda} \\ -\boldsymbol{\beta}\_{\Lambda} & \boldsymbol{\alpha}\_{\Lambda} & 1 & z\_{\Lambda} \\ 0 & 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} \mathbf{x} + \mathbf{x}\_{\Lambda} + \boldsymbol{\beta}\_{\Lambda} \cdot \mathbf{z} - \boldsymbol{\gamma}\_{\Lambda} \cdot \mathbf{z} \\ \boldsymbol{y} + \boldsymbol{y}\_{\Lambda} - \boldsymbol{\alpha}\_{\Lambda} \cdot \mathbf{z} + \boldsymbol{\gamma}\_{\Lambda} \cdot \mathbf{x} \\ \boldsymbol{z} + \boldsymbol{z}\_{\Lambda} + \boldsymbol{\alpha}\_{\Lambda} \cdot \mathbf{y} - \boldsymbol{\beta}\_{\Lambda} \cdot \mathbf{x} \\ 1 \end{bmatrix},\tag{12}$$

Then, the assembly accuracy would be

$$T\_j^{i - \omega\_\Lambda} = f\_j^i(\mathbf{K}\_i^R, \mathbf{K}\_i^A \boldsymbol{\omega}\_\Lambda^H) = f\_j^i(\mathbf{P}\_i^R, \mathbf{G}\_i^R, \mathbf{P}\_i^A \boldsymbol{\omega}\_{\Lambda'}^H \mathbf{G}\_i^A),\tag{13}$$

On this pose, if the assembly accuracy requirements are still satisfied as

$$\forall f \big( \mathsf{K}\_{i}^{R}, \mathsf{K}\_{i}^{A} \mathsf{o}\_{\Delta}^{H} \big) \in \left[ T\_{j}^{i-\min}, T\_{j}^{i-\max} \right], i \in [1, n], j \in [1, j\_{i}], \tag{14}$$

the pose is still a valid pose. The coordination space model can, hence, be expressed as

$$\mathcal{Q}\_{\mathbb{C}S} = \left\{ \boldsymbol{\omega}\_{\Delta} \middle| \forall f\_j^i \Big( \mathbb{K}\_i^R, \mathbb{K}\_i^A \boldsymbol{\omega}\_{\Delta}^H \big) \in \left[ T\_j^{i-\min}, T\_j^{i-\max} \right], i \in [1, n], j \in [1, j\_i] \right\} \tag{15}$$

where ∅*CS* is the coordination space, which is the whole pose variation space under the condition of assembly accuracy requirements.

#### **3. Assemblability Analysis Based on Coordination Space Model**

Assemblability refers to the geometric consistency of the matching geometric structures of the two assembling parts. It can be judged whether the assembly can directly be carried out by assemblability prediction.

The assemblability is good if the coordination space is greater than 0, which means at least one pose conforms to Equation (15). Otherwise, the assemblability is bad. Therefore, the assemblability analysis flow is shown in Figure 5.

**Figure 5.** Assemblability analysis process.

Firstly, the KACs are measured by a laser tracker or other digital measurement devices. Then, the coordination space model is constructed based on the assembly accuracy requirements. The volume of the coordination space is solved to judge whether it is assemblable. It will be assemblable when ∅*CS* is greater than 0. The assembly can be executed by calculating the optimal pose and aligning the parts. It will be uncoordinated when ∅*CS* is 0. Then, the assembly deviation should be analyzed and compensated to make it assemblable.

The solution process of the coordination space is based on the Monte Carlo method:


$$\mathbf{a}\mathbf{a}\_{\mathbf{d}} : (-\mathbf{x}\_{\mathbf{d}\prime} - \mathbf{y}\_{\mathbf{d}\prime} - \mathbf{z}\_{\mathbf{d}\prime} - \mathbf{a}\_{\mathbf{d}\prime} - \boldsymbol{\beta}\_{\mathbf{d}\prime} - \boldsymbol{\gamma}\_{\mathbf{d}}) \to (\mathbf{x}\_{\mathbf{d}\prime}\,\mathbf{y}\_{\mathbf{d}\prime}\,\mathbf{z}\_{\mathbf{d}\prime}\,\mathbf{a}\_{\mathbf{d}\prime}\boldsymbol{\beta}\_{\mathbf{d}\prime}\boldsymbol{\gamma}\_{\mathbf{d}\prime}),\tag{16}$$


$$\mathcal{Q}\_{\mathbb{CS}} = \frac{\mathbf{n}\_{\circ}}{\mathbf{n}\_{T}} 64 \mathbf{x}\_{\mathrm{d}} y\_{\mathrm{d}} z\_{\mathrm{d}} \alpha\_{\mathrm{d}} \beta\_{\mathrm{d}} \mathbf{y}\_{\mathrm{d}} \tag{17}$$

In the case of incoordination, the coordination space should be further analyzed. According to Equation (15), the coordination space is the intersection of KAC's constraint equations. All constraints are divided by KACs. Equation (15) will be translated to

$$\left\{ \begin{array}{c} \forall f\_{j}^{1} \big( \mathsf{K}\_{1}^{R}, \mathsf{K}\_{1}^{A} \omega\_{\Delta}^{H} \big) \in \left[ T\_{j}^{1-\min}, T\_{j}^{1-\max} \right] \big| j \in \left[ 1, j\_{1} \right], \\ \forall f\_{j}^{2} \big( \mathsf{K}\_{2}^{R}, \mathsf{K}\_{2}^{A} \omega\_{\Delta}^{H} \big) \in \left[ T\_{j}^{2-\min}, T\_{j}^{2-\max} \right] \big| j \in \left[ 1, j\_{2} \right], \\ \dots \\ \forall f\_{j}^{n} \big( \mathsf{K}\_{n}^{R}, \mathsf{K}\_{n}^{A} \omega\_{\Delta}^{H} \big) \in \left[ T\_{j}^{n-\min}, T\_{j}^{n-\max} \right] \big| j \in \left[ 1, j\_{n} \right] \end{array} \right\} \tag{18}$$

Let

$$\mathcal{O}\_{i}^{KAC} = \left\{ \boldsymbol{\omega}\_{\Lambda}^{i-KAC} \, \middle| \, \, \forall f\_{j}^{i} (\mathbf{K}\_{i}^{R}, \mathbf{K}\_{i}^{A} \boldsymbol{\omega}\_{\Lambda}^{H}) \in \left[ T\_{j}^{i-\min}, T\_{j}^{i-\max} \right] \, \middle| \, j \in [1, j\_{i}] \right\} \tag{19}$$

where ∅*KAC <sup>i</sup>* is the *KAC* coordination space formed by the assembly accuracy requirements of a *KAC*, and ω*i*−*KAC* <sup>Δ</sup> is an SDT in the *KAC* coordination space. The <sup>∅</sup>*CS* would be

$$\mathcal{Q}\_{\text{CS}} = \bigcap\_{i=1}^{n} \mathcal{Q}\_{i}^{\text{KAC}},\tag{20}$$

The relationship between the KAC coordination space and the assembly coordination space is shown in Figure 6a.

**Figure 6.** The relationship between the part and feature PFS.

Figure 6b shows the status of the *KAC* coordination space when the assembly is uncoordinated. Each area of the same color represents a *KAC* coordination space. The divided zone is named the coordination zone. The accuracy compensation method is needed to improve the assemblability.

According to Figure 6b, set the union of *KAC* coordination space as a union coordination space. It is formulated as

$$\mathcal{Q}\_{\rm UCS} = \bigcup\_{i=1}^{n} \mathcal{Q}\_{i}^{\rm KAC},\tag{21}$$

where ∅*UCS* is the union coordination space. In the union coordination space, all poses are valid for some KACs but not valid for all. Some divided zones are valid for more KACs than others, e.g., the two zones marked with 3 are better than those marked with 1 or 2. The marked number is named the coordination zone index, which is the valid KACs number in the coordination zone. If a pose in the zone marked with 3 is selected, only one KAC needs to be compensated. In this way, an assemblability optimization method is put forward by selecting a coordination zone with larger volume and KAC number. The larger volume means a better geometric consistency, and the larger KAC number means fewer KACs need to be compensated.

#### **4. Assemblability Optimization Based on the Union Coordination Space**

The accuracy compensation process is time- and effort- consuming [22] when the assemblability is poor. For example, it needs programming, clamping, tool setting, machining, loosen clamping, and other steps when finishing a KAC with cutting. Therefore, reducing the number of KACs to

be processed is an effective means to improve the assembly efficiency in many cases. The optimal pose is usually obtained under the condition of optimal assembly accuracy. If each unqualified KAC is compensated one by one under the optimal pose, more work may be needed and the assembly quality might not be good, due to the unknown assembly quality after accuracy compensation. If the assembly quality is bad after compensation, there are no alternative compensation schemes based on the least-squares method. Therefore, the assemblability optimization method is proposed to solve the incoordination problem. The key to optimize the assemblability is whether there is one or more coordination zones that can satisfy assembly accuracy requirements with fewer KACs to be compensated and a better or approximate volume of coordination space.

The coordination zone index shows the valid KACs in the certain coordination zone. The total number of all KACs is *nKAC*. The incoordination zone index shows the number of uncoordinated KACs in the coordination zone. Their relationship is

$$
\mathfrak{n}\_{\text{IZI}} = \mathfrak{n}\_{\text{KAC}} - \mathfrak{n}\_{\text{CZI}}.\tag{22}
$$

where *nIZI* is the incoordination zone index and *nCZI* is the coordination zone index. If the accuracy of uncoordinated KACs is compensated well in the coordination zone, this coordination zone will change to the assembly coordination space, as Figure 7 shows.

**Figure 7.** The coordination zone variation process by accuracy compensation.

In this way, each coordination zone can be analyzed to check whether it is good to be compensated or not. Two indicators of the coordination zone should be analyzed, one is the incoordination zone index, and the other is the volume of the coordination zone. The Monte Carlo method of Section 3 is improved to judge the state of each coordination zone one by one, and the optimal assemblability optimization schemes of the coordination zone are selected for recording.

The solution process based on the Monte Carlo method is as follows:


$$\left\{\Gamma \middle| \Gamma\_{i} = \left(\mathfrak{n}\_{\text{IZI}}, V\_{\text{CZ}}, \mathfrak{b}\_{f}, \mathfrak{s}\_{\text{\textquotedblleft}}\right), i < \Gamma\_{\text{num}}\right\},\tag{23}$$

where Γ*<sup>i</sup>* is the *i*th scheme, *nIZI* is the incoordination zone index, *V*CZ is the space amount of the coordination zone, *bf* is the information of uncoordinated KACs, *s*<sup>ω</sup> is the SDT set, and Γ*num* is the max number of the schemes.

7. Calculate the center SDT of the SDTs in the selected scheme. The assembly deviation of target features under the SDT is analyzed and the accuracy compensation is carried out.

$$
\omega\_{\Delta}^{\varepsilon} = \frac{1}{n\_s} \sum\_{i=1}^{n\_s} \omega\_{\Delta i} \tag{24}
$$

where ω*<sup>c</sup>* <sup>Δ</sup> is the center SDT, *ns* is the SDT number of *s*ω, and ωΔ<sup>i</sup> is an SDT of *s*ω. All assembly accuracies on ω*<sup>c</sup>* <sup>Δ</sup> are calculated. Then, the deviations on the excessive KACs will be compensated.

**Figure 8.** Random sampling by Monte Carlo method.

According to Equations (13) and (24), the compensation amount would be

$$C\_j^i = T\_j^{i - \omega \kappa\_\Lambda^c} - T\_j^{i - \text{cpt}} = f\_j^i(\mathbf{K}\_i^R, \mathbf{K}\_i^A \boldsymbol{\omega}\_\Lambda^c) - T\_j^{i - \text{cpt}},\tag{25}$$

where *Ci <sup>j</sup>* is the compensation amount of the *j*th assembly accuracy requirement of the *i*th *KAC*, *T <sup>i</sup>*−ω*<sup>c</sup>* Δ *j* is the assembly accuracy on the SDT ω*<sup>c</sup>* <sup>Δ</sup>, and *<sup>T</sup>i*−opt *<sup>j</sup>* is the optimal value of the assembly accuracy.

#### **5. Case Study**

#### *5.1. Space Manipulator Assembly*

The space manipulator is fixed on the spacecraft, which needs a high assembly accuracy to guarantee the stability when the spacecraft is flying. The assembly is executed by shaft and hole connectors, which are shown in Figure 9a. The connector is shown in Figure 9b. The manipulator is the align part and the spacecraft is the reference part.

**Figure 9.** Space manipulator assembly.

The upper connector is fixed on the manipulator, and the bottom connector is fixed on the spacecraft. The KACs are the assembly of the connectors. Due to the slight deformation of the spacecraft and the installation error of the bottom connectors, it is difficult for the connectors to accurately assemble at one time during the assembly of the spacecraft and the manipulator. In the original assembly process, it is necessary to try the assembly first, measure the assembly deviation of the clearance and coaxiality of each pair of connectors, make the accuracy compensation, and retry the assembly to ensure the assembly quality. The assembly takes a long time and the connectors are not convenient to be operated on the spacecraft. Therefore, the laser tracker is used to measure the connectors between the spacecraft and the manipulator. The methods in Sections 2 and 3 are taken to evaluate the assemblability based on the measurement data. The method in Section 4 is used to find the key connectors to make the accuracy compensation. The assembly is carried out after the accuracy compensation. In this way, the assembly quality is better guaranteed and the assembly efficiency is improved. The flow of the proposed method and the comparison with the original method are shown in Figure 10.

**Figure 10.** The comparison of the original and proposed assembly processes.

As shown in Figure 10, the assembly process is developed toward digital measurement and analysis. The results obtained in the actual assembly and inspection are replaced by the analysis of the measurement data. Therefore, some unnecessary assembly processes are eliminated and the possibility of repeated trial assembly is greatly reduced.

The assembly accuracy requirements of the connector are coaxiality *dr* and clearance *dc* on the matching surface, as shown in Figure 9b. The coaxiality requirement is 0.2 mm, and the clearance requirement is 0.1 mm. Assembly accuracy is compensated by gasket compensation, finishing, or position movement according to the deviation.

#### *5.2. Coordination Space Model*

The measurement of the connector is based on the measurement auxiliary tool, which is shown in Figure 11.

**Figure 11.** Measurement auxiliary tool.

After inserting the shaft into the corresponding hole, measure the four holes of the measurement auxiliary tool. The measurement data are processed as the position *<sup>p</sup>* and orientation *po*.

$$\begin{cases} \begin{array}{l} p\overrightarrow{o} = \frac{(c\_3 - c\_1) \times (c\_4 - c\_2)}{|(c\_3 - c\_1) \times (c\_4 - c\_2)|}\\ p = \frac{1}{4} \sum c\_i - l \cdot p o \end{array} \end{cases} \tag{26}$$

where *c*1, *c*2, *c*3, and *c*<sup>4</sup> are the points measured by the laser tracker.

As shown in Figure 12, the clearance *dc* and the coaxiality *dr* of a connector are

$$\begin{cases} \begin{array}{c} dr = \left| \overrightarrow{p\_2} \overrightarrow{p\_1} \cdot \sin \theta\_2 \right| \\ d\varepsilon = \left| \overrightarrow{p\_1} \overrightarrow{p\_2} \right| \cdot \cos \theta\_2 + r \cdot \sin(\theta\_2 - \theta\_1) \end{array} \tag{27}$$

where <sup>θ</sup><sup>1</sup> and <sup>θ</sup><sup>2</sup> are the angles between *<sup>p</sup>*1*p*<sup>2</sup> and *<sup>p</sup>*1*o*<sup>1</sup> or *<sup>p</sup>*1*o*2; <sup>θ</sup><sup>1</sup> can be calculated by <sup>θ</sup><sup>1</sup> <sup>=</sup> *arccos <sup>p</sup>*1*p*2· *<sup>p</sup>*1*o*1/ *<sup>p</sup>*1*p*<sup>2</sup> *<sup>p</sup>*1*o*<sup>1</sup> and, similarly, θ<sup>2</sup> can be calculated by the same way; and *r* is the radius of the matching surface, which is 15 mm. There are 20 connectors to be guaranteed at the same time.

**Figure 12.** Assembly geometric constraints analysis.

Therefore, the coordination space model is

$$\left\{ \omega\_{\Lambda} \left| \begin{array}{c} \stackrel{i\rightarrow\ \rightarrow}{\left| p\_{2}^{i-\omega\_{\Lambda}} p\_{1}^{i} a\_{\Lambda}^{j\_{1}} \sin \theta\_{2}^{i-\omega\_{\Lambda}} \right| < 0.2, \\ 0 < \left| p\_{1}^{i} p\_{2}^{i-\omega\_{\Lambda}} \right| \stackrel{i\rightarrow\omega\_{\Lambda}}{\left| \cos \theta\_{2}^{i-\omega\_{\Lambda}} + r \sin \left( \theta\_{2}^{i-\omega\_{\Lambda}} - \theta\_{1}^{i-\omega\_{\Lambda}} \right) < 0.1, \right. \end{array} \right. \right\},\tag{28}$$

where ω<sup>Δ</sup> is the random SDT based on the optimal pose derived from the least-squares method, and *pi*−ω<sup>Δ</sup> <sup>2</sup> is the parameters of Equation (27) changed by ω<sup>Δ</sup> according to Equation (12), which are listed in Equation (29):

$$\begin{cases} \begin{array}{c} \mathbf{p}\_{2}^{j-\omega\_{\Delta}} = \mathbf{p}\_{2}^{j}\mathbf{w}\_{\Delta}^{H} \\ \mathbf{p}\_{2}^{j-\omega\_{\Delta}}\mathbf{o}\_{2}^{j-\omega\_{\Delta}} = \frac{(\mathbf{c}\_{1}^{j}-\mathbf{c}\_{1}^{j})\times(\mathbf{c}\_{4}^{j}-\mathbf{c}\_{2}^{j})}{|(\mathbf{c}\_{3}^{j}-\mathbf{c}\_{1}^{j})\times(\mathbf{c}\_{4}^{j}-\mathbf{c}\_{2}^{j})|} \\ \mathbf{o}\_{1}^{j-\omega\_{\Delta}} = \arccos\left(\mathbf{p}\_{1}^{j}\mathbf{p}\_{2}^{l-\omega\_{\Delta}}\cdot\mathbf{p}\_{1}^{l}\mathbf{o}\_{1}^{l}/\mathbf{p}\_{1}^{l}\mathbf{p}\_{2}^{l-\omega\_{\Delta}}\right) \begin{vmatrix} \mathbf{p}\_{1}^{j}\mathbf{o}\_{1}^{l} \\ \mathbf{p}\_{1}^{j}\mathbf{p}\_{2}^{l-\omega\_{\Delta}}\cdot\mathbf{o}\_{2}^{l-\omega\_{\Delta}}\end{vmatrix} \end{cases} \tag{29}$$

#### *5.3. Assemblability Analysis*

Part of the raw data is listed in Table 1. All the measurement data are listed in Appendix A, Table A1.


**Table 1.** Part of the raw measurement data.

The least-squares method is taken to calculate the optimal pose and the deviations on the optimal pose. The deviations of the connectors are listed in Table 2 calculated by Equation (27).

**Table 2.** Assembly deviation prediction by least-squares method.


It can be seen that 10 connectors need to be adjusted or repaired based on the least-squares method. The coordination space is 0 at the optimal pose based on the method in Section 3, which means it cannot be assembled directly. Therefore, the assemblability should be optimized.

#### *5.4. Assemblability Optimization*

The proposed method in Section 4 is taken to find the accuracy compensation schemes. The results are shown in Figure 13.

**Figure 13.** The accuracy compensation schemes.

Figure 13 shows the accuracy compensation schemes. The first point on the X axis is the KAC quantity to be compensated. The second point is the volume of the coordination zone of the scheme. The latter ones are the number of KACs. The Y axis is the scheme number. The Z axis is the value of the X axis. Seven connectors need to be adjusted to complete assembly in scheme 1. Finally, scheme 18, which needs nine connectors to be compensated, is taken by considering the assembly quality. The coordination space of the scheme is 56d. d is the volume of the maximum pose space divided by random times. In this case, d is 2.46 <sup>×</sup> <sup>10</sup>−<sup>14</sup> mm3rad3. The KAC number to be compensated is 5, 6, 8, 10, 12, 14, 18, 19, and 20. The center SDT of the coordination zone in scheme 18 is <sup>−</sup>0.0363 mm, 0.0210 mm, 0.0098 mm, 2.52 <sup>×</sup> <sup>10</sup>−<sup>6</sup> rad, 1.64 <sup>×</sup> <sup>10</sup>−<sup>6</sup> rad,−2.35 <sup>×</sup> <sup>10</sup>−<sup>6</sup> rad .

The deviation is calculated under the center SDT listed in Table 3.


**Table 3.** The deviations of the connectors.

After simulation accuracy compensation for the above nine connectors, which is in bold and italics in Table 3 (the proposed method), the coordination space is 101d, which is greater than 0. The average coaxiality is 0.068 mm and the average clearance is 0.014 mm. The assemblability is good and the assembly can be executed directly.

After simulation accuracy compensation for the above 10 connectors, which is in bold and italics in Table 2 (the least square method), the coordination space is 9d. The average coaxiality is 0.084 mm and the average clearance is 0.014 mm. The assemblability is good but the assembly quality on coaxiality is worse.

The result shows that the proposed method will generate a better accuracy compensation scheme with less workload and better assembly quality, which improves the assemblability.

The measurement and connector adjustment process took about 8 h during the assembly. The pose adjustment process took about 2 h. Therefore, it took about 10 h in total based on the proposed method. The original assembly process took more than 20 h because the first trial assembly and accuracy compensation process cannot realize the re-trial assembly smoothly. Three or four times the assembly are needed to guarantee the assembly quality.

#### **6. Discussion**

Compared to the previous research, the major contributions in this paper are listed as follows: (1) The concept of assemblability and coordination accuracy in the design/drawing stage are extended into the measurement-assisted assembly. (2) An assemblability analysis method based on the measurement data and the coordination space model is proposed for predicting the key assembly deviations. (3) The accuracy compensation methods based on the optimal pose might lead to more workload and worse assemblability. Therefore, an assemblability optimization method is proposed for less workload and better assembly quality. In addition, the space manipulator assembly is taken as an example. The result shows that the proposed method can optimize the assemblability with less workload and better assembly quality compared to the accuracy compensation method based on the optimal pose.

The assemblability optimization method based on accuracy compensation improves the ability to detect assembly problems in advance, which will benefit the automation assembly. Further, the coordination space model and the small displacement torsor are useful for analyzing the assemblability and optimizing the tolerances in the design/drawings phase, but the assemblability optimization method is not useful. In the implementation of the method, high-precision digital measurement equipment are needed. Measurement uncertainty will affect the reliability of the final results.

Future research include evaluating the influence of the measurement uncertainty on the coordination space model. Then, the uncertainty of pose adjustment should be taken into consideration compared to the volume of the coordination space to judge the feasibility of automatic pose adjustment.

**Author Contributions:** Conceptualization, Z.C.; methodology, Z.C.; software, Z.C.; validation, Z.C.; formal analysis, Z.C.; investigation, Z.C.; resources, Z.C.; data curation, Z.C.; writing—original draft preparation, Z.C.; writing—review and editing, F.D.; visualization, Z.C.; supervision, F.D.; project administration, F.D.; funding acquisition, F.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Basic Scientific Research, grant number JCKY2016206B009 and the Marine Power Research & Development (MPRD).

**Conflicts of Interest:** The authors declare no conflict of interest.


**Table A1.** Raw measurement data.

**Appendix A**


**Table A1.** *Cont.*

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Deep Neural Network for Automatic Image Recognition of Engineering Diagrams**

#### **Dong-Yeol Yun 1, Seung-Kwon Seo 1, Umer Zahid 2,\* and Chul-Jin Lee 1,\***


Received: 11 March 2020; Accepted: 4 June 2020; Published: 9 June 2020

**Abstract:** Piping and instrument diagrams (P&IDs) are a key component of the process industry; they contain information about the plant, including the instruments, lines, valves, and control logic. However, the complexity of these diagrams makes it difficult to extract the information automatically. In this study, we implement an object-detection method to recognize graphical symbols in P&IDs. The framework consists of three parts—region proposal, data annotation, and classification. Sequential image processing is applied as the region proposal step for P&IDs. After getting the proposed regions, the unsupervised learning methods, k-means, and deep adaptive clustering are implemented to decompose the detected dummy symbols and assign negative classes for them. By training a convolutional network, it becomes possible to classify the proposed regions and extract the symbolic information. The results indicate that the proposed framework delivers a superior symbol-recognition performance through dummy detection.

**Keywords:** convolutional neural network; object detection; piping and instrument diagram; unsupervised learning

#### **1. Introduction**

Engineering diagrams (EDs) are schematic drawings describing process flow, circuit construction, and engineering device information. Among the many types of EDs, piping and instrument diagrams (P&IDs) are broadly used in the production plant industry because they contain key information of the plant, including piping, valves, instruments, control logic, and annotations. Moreover, extracting this information—for example, where or of what type the objects are—is the first step in estimating the number of elements and managing the project during its operational period. Most of the plant industries, such as oil and gas production plants, have employed large teams of engineers to manually count these entities and digitalize the information into their internal systems because there is no module available to automatically extract such information from the diagrams. For decades, these tasks have been considered inefficient and time-consuming tasks. Consequently, the demand for a module enabling an automatic engineering diagram digitalization has increased as such procedures can improve productivity and gain a competitive edge for the company in the global market.

However, there are obstacles to be overcome before applying this technology in real-world scenarios. Firstly, the symbols of P&IDs come in a diverse range of forms, with approximately 100 different types for each entity [1]. Furthermore, there is an inter-similarity between these symbols themselves. This requires the person interpreting and counting the symbolic entities to know the P&ID symbols and legend sheets. In some diagrams, it is difficult to identify—through image processing alone—a target symbol without confusing it with another symbol because there are so

many objects. Moreover, in diagrams, text information such as notes, as well as line number or size information, is presented near the symbols. This information is also often also written across the symbol; therefore, it can present another obstacle in the effective recognition of diagrams. These key challenges need to be overcome to enhance the capabilities of P&ID digitalization procedures in real-world scenarios; however, there have been few studies applied to develop an object-detection algorithm to overcome these limitations and present suitable applicability.

Within the field of machine vision for EDs, there have been a few studies that have sought to extract specific information from the diagrams. In [2], the authors present new trends on machine vision to extract various information from EDs, such as binarization, contextualization, segmentation, and recognition. One of the most popular preprocessing methods is binarization; by adopting a threshold value, this method converts an image into a binary representation, thereby removing noise and improving entity identification in the diagram. There are several methods of applying such image binarizations, including global thresholding [3], local thresholding, and adaptive thresholding [4]. For line detection in the diagrams, canny edge detection [5], Hough transformations [6], and morphological dilations have been discussed in the literature. Probabilistic Hough transform (PHT) [7], which uses a random sampling of the edge points to detect lines in images, has been applied for the robust detection of lines in engineering diagrams [8]. A shape-detection procedure, employing a consistency attributed graph (CAG) with a sliding window, was used by [9] to construct a symbol-detection procedure. As a comparable method, a recursive model of the morphological opening was implemented by [10] to identify symbols by the empty fraction of their area. For text/graphics segmentation (TGS), a connected component (CC) analysis [11] was used with size constraints in an engineering diagram [12]. In one study [13], a procedure to realize pixel-wise classification into text, graphics, and background was performed using filter banks and estimations of the descriptor sparseness. To find a specific shape in the diagrams, template matching was also used with a symbol shape incorporated as part of the prior information [14]. In [15], a threshold-based object-detection algorithm was proposed for binary images.

For graphical symbol recognition using a machine learning approach for engineering diagrams, there is very little previous research. In the past, the template matching method [16] has been famously used to find a specific shape within an image by sliding the template across the entire window. In all slides, the template calculated the similarity using convolutional estimation. However, this method has an inherent disadvantage; when it estimates the similarity within a specific region, it uses the Euclidian distance. This metric is intuitive, and the method is convenient to apply, but it cannot consider the images with large numbers of dimensions. It judges similarity using only quantitative calculations. Consequently, when it sets a high threshold value, it misses the shapes in the image owing to the stringent criteria. In contrast, when it sets a low threshold value, it detects unsuitable shapes because of its naïve criteria. This makes it difficult to find a suitable threshold value and improve its detection performance. Furthermore, almost all images in the real world contain noise or resolutions inadequate for the implementation of machine vision methods. In terms of industrial scenarios, the template matching-based method is not suitable for object detection in engineering diagrams.

In electrical EDs, circuit symbols can be recognized with morphological operations and geometric analyses [17]. A convolutional neural network can be used to recognize the symbols in hand-sketched engineering diagrams and convert them to a computer-aided design (CAD) program [18]. Adapting the Hopfield model, an iterative neural network model was implemented for symbol recognition by employing a prototype [19]. One symbol-classification method applied to P&IDs considered class decomposition using k-means clustering [20]. Fully convolutional networks (FCN), which are an end-to-end network for pixel-wise prediction, were first applied as object detection for P&IDs in [21]. In [22], the author proposed a method to extract various objects, including symbol, characters, lines, and tables in a P&ID, using a machine vision method containing deep learning architecture. To reduce human effort while validating the CAD documents, the P&IDs attributed by the graph form were trained by a neural network and predicted the components vector which represents the diagram flow

in [23]. However, these previous approaches toward symbol recognition for P&IDs could not control the detection quality effectively for the symbols in the diagrams. Their applicability to real diagrams could not be confirmed because unpredictable objects were detected during the region proposal. This makes it difficult to recognize the target symbol owing to the characteristics of engineering diagrams made up of quasi-binary components.

In object detection, certain algorithms have delivered remarkable performances in recent years. The appearance of convolution networks has heralded significant improvements in image-classification problems. Even so, unlike basic image classification, object detection requires the solution not only of multi-labeled classification problems but also the bounding boxes for proposed regions in a digital image. To achieve this, networks for object detection employ a region proposal network (RPN), which plays the role of finding the symbol-region candidates, before classifying images. In a region-based convolutional neural network (R-CNN) [24], a selective search algorithm [25] was used for the RPN, and the boxed images were fed to the network for classification. RPN extracts numerous boxes from an image, considering the colors, scale, boundary, etc., of the object. The proposed regions are reshaped before being fed into the convolutional network for image classification. However, R-CNN has inherent limitations; it is expensive and slow. All processes in R-CNNs, from the RPN to the convolutional network, render the model inefficient and slow to detect objects in an image in real-time. Advanced models such as Fast R-CNN [26] and Faster R-CNN [27] improve network performance and speed by including the RPN in the neural network. In fast R-CNN, by employing a selective search algorithm as the RPN and implementing it into the neural network, it becomes possible to combine the different procedures into one end-to-end network. In Faster R-CNN, to reduce time consumption in the selective search algorithm, the algorithm was replaced by a combined neural network, which made the detection much faster. Many other state-of-the-art algorithms are being proposed as next-generation object identification strategies [28,29]. However, these have focused on problems such as the improvement of model accuracies for colored images, improvement of detection speeds, and image segmentation (i.e., the process of partitioning the image into a set of pixels with multiple segments). However, engineering diagrams have a characteristic difference from colored images; they are an almost-binary component matrix with a specific size and shape for each symbol. This makes it difficult for a model to classify symbols using only limited information.

As an advanced application of object detection in EDs, we propose an R-CNN architecture containing clustering. For the RPN, we implement a sequential image-processing method that is modified for two types of target symbols: valves and instruments. By modifying the image-processing method for region detection, we propose candidate symbol regions using size-based detection. After the RPN, we get the symbols, but we also get the meaningless regions inevitably, such as truncated line, curve, noise, and so on, we call these 'dummy'. The detected dummy images are decomposed by unsupervised learning methods, and negative classes are assigned to them for image classification. For images containing positive classes, which are our target symbols, the dataset is augmented with padding-block. Through a simple convolutional network, the multi-class classification model is trained and applied to new diagrams for the model test.

In this research, we propose a model based on an R-CNN architecture that features dummy image clustering. A sequential image-processing method is used for the RPN, instead of a selective search algorithm. After the RPN, the dataset is constructed by positive classes, and dummy clustering is applied to treat the unwelcomed detections as negative classes; thereby improving the classification performance of the convolutional network. The remainder of this paper is organized as follows. Section 2 provides the methodology for the extraction of target symbols from P&IDs. Based on our proposed method, the region proposal and symbol-recognition results are discussed in Section 3. Finally, we conclude the paper in Section 4.

#### **2. Materials and Methods**

Here, we propose our R-CNN framework for recognizing graphical symbols in P&IDs, as shown in Figure 1. There are two main types of graphical symbols targeted in this study: valves and instruments. These symbols have characteristics such as size and shape, as shown in Figure 2. Using these characteristics, we construct an RPN by modifying several image-processing techniques. There are two types of proposed regions: symbol and dummy. For data annotation of the symbols representing positive samples, a P&ID symbol and legend sheet are referred to, which provide the standard set of shapes and symbols for documenting the diagram.

**Figure 1.** A framework summary for symbol recognition of the piping and instrument diagram (P&ID).


**Figure 2.** Graphical characteristics of the target symbols.

For the dummy images representing negative samples, two unsupervised learning algorithms: k-means clustering and deep adaptive clustering (DAC), are used to analyze their hidden patterns and assign classes. After the annotation is completed, data augmentation is applied to generate additional information about the symbols. A convolutional neural network (CNN) is used to classify the symbols in this research owing to its superiority in local feature extraction [30]. After training the network, we apply it to another diagram in the same project and verify the results.

#### *2.1. Data Sets*

To implement the proposed framework and validate its performance, we use 10 pages of P&IDs from a real project. The resolution of the diagrams is 300 dpi in A3 size; thus, they contain approximately 4000 × 3000 cubic pixels. Of the 10 pages, we take seven and apply the region proposal method to construct our dataset of the proposed regions, which contains both positive and negative samples. After augmentation by padding (100 × 100 pixels), this is fed into a simple CNN.

We construct and compare three models based on the type of data that they use—positive samples only (P), positive with negative samples through k-means clustering (PN\_Kmeans), and positive with negative samples through deep adaptive clustering (PN\_DAC). To investigate and test our models, they are coded using Python with a Tensorflow backend. We also maintain the same computational conditions using NVIDIA TITAN V with 8 GB GDDR5.

#### *2.2. Region Proposal*

Region proposal is the process to extract candidate regions of symbols. Instead of apply the selective search algorithm, we build a customized process to extract the candidate regions in EDs as given in Figure 3.

**Figure 3.** Procedure followed by the proposed region proposal network (RPN).

There are several points of incongruence in the selective search algorithm. For the proposal of candidate regions, the selective search algorithm uses the traits of an image, such as its color and boundary. However, this is not appropriate for the detection of an object in a binary image, such as an engineering diagram. To investigate the most-suitable algorithm for an engineering diagram, we implement sequential image processing to create proposal regions; this detects the potential target symbols by their characteristics. The sizes and shapes of plant symbols are specified in the diagram; therefore, it is possible to modify the image processing technique for each type of target symbol—valves or instruments—as they have a set size and aspect ratio in copies with identical resolutions. Thus, using copies of the input image, the characteristics of each target symbol can be reflected and sought for in each step. To modify the progress, we divide it into four parts: (1) image binarization, (2) non-target removal, (3) morphological transformation, and (4) CC analysis, as shown in Figure 3.

In image binarization, the adaptive threshold method [31] is used to reduce the noise present in the input images and convert them into a binary representation. Comparing a pixel against the average of those surrounding it preserves hard contrast lines and discards soft gradient changes.

In the non-target removal step, we remove the obstacles for the detection of each symbol. In the case of valve detection, other symbols such as lines, instruments, and pipe fittings are considered obstacles to clear detection. This process is employed as an intermediate stage to remove the obstacles and reduce the number of meaningless detections in the region-proposal step. For line removal, dilation kernels are used in the horizontal and vertical directions, with the structures being (1 × p) and (q × 1), respectively. The kernel parameters p and q are adjusted by considering the size of the symbol. In this study, we use a length similar to the shortest side of the symbol. For non-target symbol removal, a CC analysis algorithm [11] and Hough circle algorithm [32] are used to find the contours of the non-targeted objects such as instruments and pipe fittings.

In the morphological transformation step, "Closing," which is derived from the basic operations of erosion and dilation, is commonly used to enhance object outlines and small cover-up holes in the image [33]. Through this closing method, the floating objects retained from previous steps protect those background regions that have a similar shape to their kernel, while deleting all other background pixels [34].

Finally, to propose regions for the candidates of the target symbol, CC analysis is used with the constraints of symbol size and aspect ratio. The algorithm analyzes the topological structure of binary images. At the level of individual pixels, it considers 4-(8-) neighboring regions for the connected cases. We assume each symbol size and aspect ratio as prior knowledge in the detection. Figure 2 presents the schematic procedure of our region-proposal method. The image processing for region proposal plays a role in reducing the total number of proposed regions, and it adjusts the detection of undesirable objects called dummies.

After the region proposal step, we construct a dataset for the classification network as given in Figure 4. There are two types of proposed regions: symbol and dummy. Symbols are our positive samples; they are the gate valves, check valves, sensors, etc. Dummy entities have unpredictable shapes and sizes, are not within our interest, and make it difficult to classify symbols through the machine learning algorithm. Therefore, they are considered to be negative samples. For the positive dataset, we use the P&ID symbol and legend sheets, which provide a standard set of symbol shapes and legends for documenting diagrams and assign the symbol images for each class, such as gates, balls, globes, checks, etc. From the proposed regions, the positive samples are manually classified into 10 classes according to their shape and function.

**Figure 4.** Data annotation: positive/negative samples.

#### *2.3. Dummy Clustering*

Aside from the positive samples, numerous images remain from the proposed regions, which are called dummies. They consist of curved lines, arrow shapes, revision clouds, or cut entities, as shown in Figure 3. It is difficult for the region proposal method to control the detection of dummy entities because the diagrams are quasi-binary representations and consist of many entangled lines and entities. In this research, we assign the negative classes of the detected dummies to improve the classification performance and consider the applicability of the procedure to real projects.

Dummies are obstacles for the classification model seeking to identify target symbols from a pool of proposed regions. P&ID is a type of grayscale image that is composed of only one channel, therefore, there is an arbitrary limitation to classifying the proposed regions when only using the positive samples from the target data. Furthermore, in the region-proposal network, the patterns of

the detected dummies are unpredictable because it is difficult to erase all non-target entities during image processing. These patterns, such as shape, size, and detection frequency, are uncertain in every diagram, and this makes it difficult for the model to identify the target from the pool of detection images. Therefore, we assign negative samples to the classification models with unsupervised learning algorithms. To decompose the pool of dummy images and assign the class as a negative sample, we apply two unsupervised learning algorithms—k-means clustering and deep adaptive clustering (DAC). K-means clustering is a basic unsupervised learning algorithm. It is an iterative method to locating k-centroids in the dataset [35]; it locates them by optimizing the position of each centroid, based on the L2 norm in the feature space, as shown in Equation (1):

$$\operatorname{argmin}\_{\mathbb{C}} \sum\_{i=1}^{K} \sum\_{\mathbf{x}\_{j} \in \mathcal{C}\_{i}} \|\mathbf{x}\_{j} - \mathbf{c}\_{i}\|^{2} \tag{1}$$

$$X = \mathbb{C}\_1 \cup \mathbb{C}\_2 \cdots \cup \mathbb{C}\_{K'} \; \mathbb{C}\_i \cap \mathbb{C}\_j = \phi \tag{2}$$

The quantity *x* represents a pool of unlabeled data, and *ci* is a centroid of the *i*-th cluster, *Ci*.

DAC is also applied to decompose the hidden patterns of the dummies with an advanced method [36]. It is one of the state-of-the-art algorithms for the image-clustering problem that uses a convolutional architecture and cosine distance to measure the similarity of pairwise images with adaptive parameters. It delivers superior performance in image clustering owing to its adaptive-learning algorithm. The network solves the image-clustering problem as a binary pairwise-classification problem. The flowchart is presented in Figure 5.

**Figure 5.** Flowchart of the deep adaptive clustering (DAC) algorithm [36].

Initially, unlabeled images are input to the convolutional network to generate a provisional latent vector for the image. Using the latent features, cosine similarities between pairwise images *xi* and *xj* are calculated; then, a confusion matrix is constructed for every batch. The network objective function is defined by:

$$Min\_0 \mathcal{E}(\boldsymbol{\theta}) = \sum\_{i,j} L(r\_{i\mid \prime}, l\_i \cdot l\_j)\_{\prime} \tag{3}$$

$$\left\|\mathbf{s}.t.\,\forall i,\,\left\|l\_{i}\right\|^{2} = 1,\,\text{and}\,l\_{i\text{l}} \ge 0,\,h = 1,\,\dots,\,k\right.\tag{4}$$

where *rij* is the unknown binary variable; if the pair of input images are in the same cluster, then *rij* = 1, otherwise, *rij* <sup>=</sup> 0. · <sup>2</sup> represents the L2 norm of a vector, and *lih* represents the *<sup>h</sup>*-th element of the label feature of the k-dimensional latent vector *li*. As the cosine similarity of the input image pair

can be formulated by *li*·*lj*, the objective function of DAC is expressed by the loss between *rij* and *li*·*lj*. The expression for *L rij*, *li*·*lj* is formulated as follows:

$$L\left(r\_{i\bar{j}\cdot}, l\_{i\cdot}l\_{\bar{j}}\right) = -r\_{i\bar{j}} \log\left(l\_{i\cdot}l\_{\bar{j}}\right) - \left(1 - r\_{i\bar{j}}\right) \log\left(1 - l\_{i\cdot}l\_{\bar{j}}\right) \tag{5}$$

However, the unknown variable *rij* is prior information. Thus, an adaptive parameter λ is applied for the stepwise threshold value; we use μ(λ) and *l*(λ) as the values for selecting similar (or dissimilar) image pairs.

$$r\_{ij} = \begin{cases} 1, \text{ if } l\_i \cdot l\_j \ge \mu(\lambda) \\ 0, \text{ if } l\_i \cdot l\_j \le l(\lambda), \text{ i.e.} \\ \text{None, otherwise.} \end{cases} \tag{6}$$

In the clustering process, the value of λ starts at a specific value and gradually increases. Besides this, the relationships μ(λ)∝λ, *l*(λ)∝λ, and *l*(λ) ≤ μ(λ) are set in the algorithm. After finishing a batch process, the parameter λ is also updated by the gradient descent algorithm.

$$Min\_{\lambda}E(\lambda) = \mu(\lambda) - l(\lambda) \tag{7}$$

$$
\lambda := \lambda - \eta \cdot \frac{\partial E(\lambda)}{\partial \lambda} \tag{8}
$$

Here, η represents the learning rate of λ. Using this adaptive modification of the parameter λ, the algorithm performs a stepwise selection between the pair images with increasing λ. The performance of the DAC is detailed for various datasets in [36]; it delivers the best performance in a binary image clustering problem, such as MNIST when compared against other clustering methods. Therefore, as an advanced method to decompose dummy images, DAC is applied in this research, and the results are compared with those of the k-means clustering.

The detailed architecture of DAC is summarized as follows. After the input images are padded by 100 × 100, we use six convolutional layers with a (3 × 3) kernel size, (1 × 1) stride, ReLU (Rectified Linear Unit) activation, and padding of the same structure in this network. The number of filters in each layer are 64, 64, 64, 128, 128, and 128, respectively. A max-pooling operation is applied with a (2 × 2) kernel and (2 × 2) stride. In fully connected layers, hidden units contain 128 and 64 nodes with ReLU activation. In all the layers, batch normalization is used to prevent the outputs of the hidden nodes from fluctuating. For adaptive learning, we set the selection-control equations according to Equations (9) and (10):

$$
\mu(\lambda) = 0.95 - \lambda \tag{9}
$$

$$l(\lambda) = 0.455 + 0.1 \cdot \lambda \tag{10}$$

There are 451 instances of dummy images from the seven pages of P&IDs. Using these two algorithms, the hidden patterns of the dummy pool are identified; then, we automatically assign them into k classes as negative samples. In this research, the value of k is fixed at 13. The optimum value of k is also an issue in clustering problems; however, we only focus on the effects of assigning a negative class for classification networks.

Table 1 presents the results of the data structure with positive and negative classes. Based on the clustering results, 13 negative classes are constructed, along with 10 positive classes. After annotating the data, which contains a total of 23 classes for the proposed regions, data augmentation [37] is applied to enhance the information in each class. To augment data, we apply two methods, central movement, and rotation. Central movement means that the extracted images from the RPN are padded by 100 × 100, before entering the network. We implement it moves 1 × 1.3 × 3 pixels around the center of the image to catch a located symbol a little sideways for the same class. The rotation is also applied because some valves exist in a rotated form in the diagrams, so, in this study, only 45, 90, 135, 180, 225, 270, and 315 rotation angles were implemented to catch the rotated symbol.


**Table 1.** Data annotation and augmentation.

#### *2.4. Convolutional Network*

Several machine-learning methods can be applied to the image-classification problem, including the support vector machine, random forest, and neural network-based models; however, we implement a simple convolutional neural network as our classification model to extract local information of the image data by convolutional and max-pooling filters [37].

The detailed model structure is presented in Figure 6. We construct three convolution layers and two fully connected layers in the network. The number of convolution filters in each layer is 64, 128, and 256, respectively. A kernel size of (3 × 3), a stride of (1 × 1), and a max-pooling layer with (2 × 2) filters are used for local feature extraction. Fully connected layers consist of 256 and 23 units, and the ReLUactivation function is used throughout our model, except for the end unit, wherein a softmax function is used. For generalization of the model, the dropout method [38] is applied, which is set at 0.7. The purpose of the dropout is to prevent overfitting problems in the neural network-based model by applying a zero forward-direction propagation value stochastically to every layer.

**Figure 6.** R-CNN scheme for symbol recognition of P&ID.

#### *2.5. Evaluation Metric*

Considering the target symbols in the diagrams and the requirements of practical applications, we suggest two metrics for validating the proposed framework-symbol recognition rate (*SR*) and dummy detection rate (*DR*), using a constant confidence threshold of 0.7.

$$\text{SSR (\%)} = \frac{\text{The Number of Correct Reynolds}}{\text{The Number of Squares in the diagram}} \times 100\tag{11}$$

$$DR\left(\%\right) = \frac{\text{The Number of Dummies Consfused with Symbols}}{\text{The Number of Model Predictions}} \times 100\tag{12}$$

*SR* is the number of correctly recognized symbols divided by the number of symbols in the diagrams. It describes to what extent the model correctly detects the target symbols in the diagrams. *DR* is calculated by dividing the number of dummy images confused with symbols by the number of model predictions. It also describes the capacity of our model to distinguish dummies from symbols. If the model is well-trained in object detection for P&IDs, the value of *SR* will be large, whereas the value of *DR* will be small. In this study, the models are validated and compared with each other using these two metrics.

#### **3. Results**

#### *3.1. Region Proposal Results*

Figure 7 summarizes the region proposal results. The target symbols-valves and instruments were well-detected in these results. We implemented a customized procedure for each target symbol and integrated the proposed regions into one diagram. All the targets in the diagram were detected using image processing. For each target, the image processing was set to modify the overlapping contours in the detected regions.

Sine the proposed regions were detected by size constraints in the contour method—that is the CC analysis—there were unwelcomed images in the resulting diagram. These had a similar size to the target and represented sliced lines, the edges of instruments, entangled lines, etc. To reduce the number of dummy detections, the size constraints were used to customize the image processing for each target by adopting the target size as prior knowledge. The main purpose of the region proposal was to identify the candidate regions where the target symbol might exist; therefore, a noteworthy advantage of the process is that we are not required to focus on making the number of candidate symbols as small as possible; they must be detected conservatively and passed into the convolutional network for target identification.

**Figure 7.** Sample results of the region proposal network (RPN).

From these proposed regions, we obtained a pool of images containing symbols and dummies. In the P model, only the symbol data are constructed as the dataset for the classification model. On the other hand, in the PN models, both symbols and dummies are incorporated into the model. To assign classes to the dummies, a series of detected dummies was decomposed through the clustering

algorithms. Consequently, the effects of negative classes on symbol recognition in engineering diagrams were analyzed; these are described in the following section.

#### *3.2. E*ff*ects of Negative Classes*

First, we only investigated the positive samples to test the performance of the model. The model recognized symbols in the test diagrams but could not distinguish dummy images from the proposed region. This demonstrates that the model, which is only trained with positive samples, can distinguish only symbols. We could observe that both PN models-k-means and DAC were superior to the P model in terms of target-symbol classification from the proposed regions. In Figure 8a, the symbols were well recognized by the P model, but dummies were also detected in the results. This means that the model, which was trained only on positive data, had a weakness in identifying negative samples as false. In contrast, PN models such as Figure 8b exhibited strong discriminative performance between symbols and dummies. The dummies confused with the check valves and gate valves were filtered out by the PN models. These results indicated that the assignment of a negative class for classification gives the model the ability to effectively identify symbols from the pool of binary component images.

**Figure 8.** Sample results from the (**a**) P\_model, (**b**) PN\_Kmeans model.

We verified this statement in Table 2. In terms of SR—the extent to which the model could recognize the symbols in the diagram—all the models delivered good performance of over 96%. The PN\_DAC model outperformed the other two, with 98.08%. This suggests that in the PN models, there was enhanced ability to classify targets through the assignment of a negative class.



DR demonstrated the remarkable ability of both PN models. In the P model, 42.3% of the dummy images from the test diagrams were confused with our target symbols. This meant that the P model was incapable of identifying what images represented a genuine symbol. Due to the characteristics of EDs, i.e., their binary representation, it was difficult for the model to recognize them. In the latent space of the convolutional network, the latent features of both the dummies and the symbols were confused with each other under the P model because it does not possess any information concerning negative images. Hence, we conclude that negative classes are required for object-detection algorithms in EDs.

Compared to the P model, the models that consider negative samples achieve a significant reduction in the dummy detection rate. Binary images contain only a limited amount of information. Althoughmost images in real-world applications consist of three channels-red, green, and blue-engineering diagrams consist of one channel-grayscale. They feature only one channel, which consists of quasi-binary components; hence, there is limited information available in the image, such as local features or pixel intensity. In this respect, we can say that for engineering diagrams to effectively recognize plant symbols and discard the detected dummies from the proposed regions, a dataset containing positive and negative classes is required. Consideration of the negative results yields the additional model information through which candidates can be assessed effectively.

#### *3.3. E*ff*ect of Clustering Methods*

In Table 3, it is shown in a confusion matrix to depict the performance of PN\_DAC model.

The PN\_DAC model exhibits the optimum performance in SR and DR, with 98.08% and 0.39%, respectively. Though both PN models-k-means and DAC-had a low dummy detection rate, the PN\_DAC model recorded a lower score in dummy detection than the PN\_Kmeans model by approximately 1%. This resulted from the differences between the image clustering methods. The k-means clustering is an iterative algorithm based upon Euclidian distance, which represents a simple quantitative distance between entities in feature space. It does not consider the direction of the feature. Consequently, the algorithm is too weak to construct with high-dimensional data such as that contained within image representations. In contrast, the PN\_DAC model obtains the latent features of the data by performing efficient feature extraction using a convolutional network. As shown in Table 3, most of the confusion is created among these symbol classes, except for the three-way valves, ball valves, and sensor symbols. We also calculate F1 scores for each symbol in PN\_DAC model, as given in Table 4. By solving the pairwise binary classification problem with adaptive parameters, the model delivered good performance that could be interpreted as a good analysis of the hidden patterns in the regions. DAC has configured the negative class to make it easier to distinguish between the classes, thereby increasing its performance by entering well-defined data into the model.


### **Table 3.** Confusion matrix: PN\_DAC model.


**Table 4.** F1 score for each class (PN\_DAC).

#### **4. Conclusions**

In this study, an R-CNN for engineering diagrams was proposed, taking negative classes into account. For an RPN, sequential image processing was modified for each target—valve and instruments. To annotate the negative class for the dummy images, two unsupervised learning algorithms–k-means and DAC–were applied to decompose the hidden patterns of the dummies, and assign negative classes. A simple convolutional network was used as the classification model because of its superior characteristics in terms of local information extraction from the images.

There were three types of datasets used for the classification problem—positive (P model), positive with negative through k-means (PN-Kmeans model), and positive with negative through DAC (PN-DAC model). Compared to the P model, both k-means and DAC had relatively low dummy detection rates of 1.35% and 0.39%, respectively, because the negative class from the unsupervised algorithm improved the model's ability to distinguish dummies from the symbols in the diagrams. Moreover, DAC was a superior algorithm for decomposing binary representations, as the PN-DAC model had superior performance, in which the symbol recognition rate (SR) and the dummy detection rate (DR) were 98.08% and 0.39 %, respectively.

From these results, we can verify that the proposed model meets the applicability and practicality criteria for P&ID object detection algorithms. The algorithm's negative sample detection reduces due to dummies, which makes its application to real projects difficult. This object detection algorithm is expected to contribute to the automatic digitalization of engineering diagrams. Regarding further work, state-of-the-art algorithms for object detection, such as Faster R-CNN, You Only Look Once(YOLO)\_v3, and Single Shot Multi-Box Detector (SSD), could be modified to suit engineering diagrams. For real-world applications, a tiny-object detector also would be useful as a plant symbol recognition model.

**Author Contributions:** Conceptualization, D.-Y.Y.; methodology, D.-Y.Y.; software, D.-Y.Y. and S.-K.S.; validation, D.-Y.Y. and S.-K.S.; formal analysis, D.-Y.Y. and S.-K.S.; investigation, D.-Y.Y.; resources, D.-Y.Y.; data curation, D.-Y.Y.; writing—original draft preparation, D.-Y.Y. and S.-K.S.; writing—review and editing, D.-Y.Y. and S.-K.S.; visualization, D.-Y.Y.; supervision, C.-J.L. and U.Z.; project administration, C.-J.L. and U.Z.; funding acquisition, C.-J.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the Chung-Ang University Research Grants in 2018 and Seoul R&BD Program (20191471).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
