# Sensitivity Analysis Based on Polynomial Chaos Expansions and Its Application in Ship Uncertainty-Based Design Optimization

Authors: Xiao Wei; Haichao Chang; Baiwei Feng; Zuyuan Liu
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

In order to truly reflect the ship performance under the influence of uncertainties, uncertainty-based design optimization (UDO) for ships that fully considers various uncertainties in the early stage of design has gradually received more and more attention. Meanwhile, it also brings high dimensionality problems, which may result in inefficient and impractical optimization. Sensitivity analysis (SA) is a feasible way to alleviate this problem, which can qualitatively or quantitatively evaluate the influence of the model input uncertainty on the model output, so that uninfluential uncertain variables can be determined for the descending dimension to achieve dimension reduction. In this paper, polynomial chaos expansions (PCE) with less computational cost are chosen to directly obtain Sobol' global sensitivity indices by its polynomial coefficients; that is, once the polynomial of the output variable is established, the analysis of the sensitivity index is only the postprocessing of polynomial coefficients. Besides, in order to further reduce the computational cost, for solving the polynomial coefficients of PCE, according to the properties of orthogonal polynomials, an improved probabilistic collocation method (IPCM) based on the linear independence principle is proposed to reduce sample points. Finally, the proposed method is applied to UDO of a bulk carrier preliminary design to ensure the robustness and reliability of the ship.


---

## Body

## 1. Introduction

Ship uncertainty-based design optimization (UDO) is a method of optimizing the design space according to the requirement of the robustness and the reliability under the influence of the uncertainty [1, 2], including robust design optimization (RDO) and reliability-based design optimization (RBDO). Papanikolaou et al. [3] considered uncertainties related to ship’s seakeeping responses and wave induced loads and presented their recent advances in modelling the combined hydrodynamic responses of ship structures using cross-spectral combination methods and using uncertainty models for the development of modern decision support systems as guidance to ship’s masters. Recently, Plessas and Papanikolaou [4] developed the procedure which may be used as a Decision Support Tool (DST) for interested ship investors including a Life Cycle Assessment (LCA) model that takes into consideration uncertainties of the most dominant economic parameters in tanker ship operation to demonstrate the optimization of alternative hull/engine/propeller setups for a defined operational scenario of a tanker that includes ship’s operational profile in calm seas and representative weather conditions, as well as all relevant safety and efficiency regulatory and technical constraints. Hou et al. [5] took Autonomous Underwater Vehicle (AUV) as the research object, the resistance as the optimization objective, and the characteristic parameters of AUV as the uncertainty and then evaluated the performance of AUV with the 6 sigma design criterion to ensure the operation reliability of AUV under the environmental interference. After that, 4 different forms of uncertainty and their applications in the RDO were considered [6], where the standard Wigley ship model was taken as the research object and EEOI of the whole route was used as the optimization objective to verify feasibility and superiority of this method. The same method was also applied in the hull form optimization of a bulk carrier [7].In addition, Matteo Diez and his team have led in the ship uncertainty-based design optimization for years. They applied the PSO algorithm to the robust optimization concept design of a bulk carrier [8] and proposed a two-stage ship design optimization method considering the uncertainty [9]. Then they systematically demonstrated RDO, introduced an integral form uncertainty quantification method, and compared various robust optimization forms to verify the superiority of the method through a robust optimization example of a bulk carrier [10]. For the epistemic uncertainty, a Bayesian method was proposed to quantify this kind of uncertainty while the probability density function was still used to model the aleatory uncertainty [11]. Subsequently, he proposed a variable-accuracy metamodel-based architecture to improve the efficiency of multidisciplinary robust optimization design and applied it to the hydrostatic resistance optimization of DTMB5415 model [12]. In order to accurately establish the probability distribution of uncertain variables, Diez and his team obtained a large amount of operational data through a two-month data collection of five sister ships and derived a probability distribution of displacement and speed [13], where some modifications of a small portion of the hull were proposed in order to increase significantly the performances of the hull and decrease the operative cost of the ship.When the ship design process is tackled, uncertain variables of the mission profile, main dimensions, etc. need to be defined. In order to truly reflect the ship performance under the influence of uncertainties, the designer gradually considers more and more uncertainties in the ship UDO process to obtain a more robust and reliable solution. It should be acknowledged that more uncertainties are helpful to simulate the real operating environment, but the high-dimensional problems that come with it also raise considerable challenges. For UDO, it is necessary to quantify the uncertainty for every case; therefore, it is obvious that adding another uncertain parameter will definitely increase the calculation burden, resulting in low optimization efficiency. Sensitivity analysis (SA), as a dimension reduction technology, is a feasible way and has been adopted to alleviate this problem. SA can qualitatively or quantitatively evaluate the influence of the model input uncertainty on the model output, so that uninfluential uncertain variables can be determined for the descending dimension to achieve dimension reduction. SA methods can be divided into local sensitivity analysis (LSA) method and global sensitivity analysis (GSA) method. The former investigates effects of variations of input factors in the vicinity of nominal values, whereas the latter aims to quantify the output uncertainty due to variations of the input factors in their entire domain. GSA does not require the model to be a linear system like LSA and it also includes investigation of interactions between model parameters, which is currently more widely applied.For GSA methods, there are Fourier amplitude sensitivity test (FAST), extended FAST (EFAST), random balance design (RBD), and Sobol' global sensitivity method based on variance decomposition, etc. Among several GSA methods, GSA with Sobol’ sensitivity indices is herein of interest [14], which belongs to the broader class of variance-based methods. These methods do not assume any kind of linearity or monotonicity of the model and rely upon the decomposition of the response variance as a sum of contributions of each input factor or combinations thereof, that is, partial variance. Various methods have been investigated for computing Sobol’ indices; among them, Monte Carlo (MC) method is still the most commonly used method to sample and estimate the integral in indices; however, a prominent drawback of this classical method is that it often requires thousands of model evaluations, which requires substantial computational resources and makes it impractical in terms of time consumption. Therefore, it is obviously infeasible to use MC method directly. In other words, a method as accurate as the MC method but with lower computational cost is desired.In recent years, PCE has gradually been introduced in this field. Using PCE method to perform SA in UDO, the coefficients of the polynomial can not only obtain the stochastic property of outputs required (mean, standard deviation, skewness, and kurtosis) in UDO, but also directly obtain Sobol' sensitivity indices. In other words, once a PCE representation is available, calculating Sobol’ indices is only the postprocessing of coefficients of PCE; therefore, indices can be calculated analytically at almost no additional computational cost. It was originally introduced by Sudret [15] and applied to replace the traditional Sobol’ MC method to calculate Sobol’ indices. In the field of the stochastic modelling of subsurface flow and mass transport, PCE has been widely used and proved to be able to perform SA comprehensively and reliably at low computational cost. Fajraoui et al. [16] and Younes et al. [17] applied a PCE-based global SA to flow and mass transport in a heterogeneous porous medium and established the transient effect of uncertain flow boundary conditions, hydraulic conductivities, and dispersivities; in the radionuclide transport simulation in aquifers, Ciriello et al. [18] analyzed the statistical moments of the peak solute concentration measured at a specific location; Sochala and Le Maitre [19] considered uncertain effects of soil parameters upon 3 different physical models of subsurface unsaturated flow; Formaggia et al. [20] applied PCE-based sensitivity indices to evaluate the influence of the uncertainty in hydrogeological variables on a basin-scale sedimentation process; Deman G. et al. [21] combined the computation of Sobol’ indices with sparse PCE method to effectively analyze the influence of input parameter uncertainty of the ground water over the life cycle on model outputs.Compared with MC method, PCE method can greatly reduce the computational cost; however, it still has the potential to further reduce the amount of calculation. In the process of solving the polynomial coefficient of PCE, probabilistic collocation method (PCM) is used instead of the common statistical methods, such as Latin hypercube sampling method, which select a large number of sample points to maintain the calculation accuracy. Generally, the number of sample points should be greater than the number of undetermined coefficients. However, different from statistical methods, due to the orthogonal property of the polynomial, input sample points of PCM are not randomly selected, but according to certain rules, which may lead to fewer sample points. Xiu et al. [22] discussed the use of PCM to solve elliptic stochastic partial differential equations; Foo et al. [23] used PCM to study the three-dimensional problem with random loads and the spatial variability material; Li et al. [24] applied PCM to the simulation of the ground water seepage and the solute transport; Isukappali [25, 26] thought that the number of collocation points should be twice as many as the number of undetermined coefficients and then used the regression method to calculate the coefficients. Unfortunately, Jiang et al. [27] found that sometimes increasing the number of collocation points to more than 8 times the number of the undetermined coefficients still cannot maintain the calculation accuracy. In order to deal with the unbalance between the number of sample points and the calculation accuracy, an improved probabilistic collocation method (IPCM) based on the linear independence principle is adopted to give the optimal number of collocation points by comparing the rank of the coefficient matrix with the number of collocation points, which reduces collocation points for solving polynomial coefficients of PCE, thereby improving the efficiency of the uncertainty-based optimization.The paper takes the ship uncertainty design as the research object, analyzes the influence of multiple random variables on the optimization target and the constraint, and proposes a ship uncertainty design based on the multidimensional polynomial chaos expansion method. It finally applies it into the practical engineering optimization design of a bulk ship, to ensure the robustness and reliability of the ship.The author has done some work about ship UDO [28]. In contrast to the previous work, this paper uses Sobol’ indices based on PCE for dimension reduction, reduces the sample points by using an improved probability collocation method, and solves failure probability via maximum entropy method (MEM) and is organized as follows. In Section 2, a brief outline of uncertainty-based design optimization and its formulation are presented. Section 3 introduces PCE theory and summarizes the process of solving the polynomial; moreover, it provides concepts of SA with Sobol’ indices and the computation of those indices using PCE. In Section 4, the basic theory of PCM is introduced and its improvement based on the linear independence principle is given with numerical examples to verify its feasibility. Finally, the proposed method is applied to the uncertainty-based optimization of the bulk carrier preliminary design in Section 5, whereas final conclusions are shown in Section 6.


## 2. Uncertainty-Based Design Optimization

According to different design requirements, uncertainty-based design optimization can be divided into robust design optimization, reliability-based design optimization, and reliability-based robust design optimization.


### 2.1. Robust Design Optimization (RDO)

Compared with the traditional design optimization, RDO considers the interference of external factors on outputs, whose solutions are not easily perturbed by external factors.As shown in Figure1, point 1 is the optimal point of deterministic optimization (DO) while point 2 is the optimal point of RDO. If the influence of uncertain factors is not considered, the objective function value of point 1 is minimum, which should be the optimal solution of objective function. However, when the design variable x perturbs in the range of ±Δx, the perturbation ΔF1 of objective function at point 1 violates constraints, resulting in failure of the product. The perturbation ΔF2 is obviously less than ΔF1, and all the objective outputs are within the constraint. Therefore, the robust optimal solution is a relatively stable optimal solution.Figure 1
Deterministic and robust optimal solutions.RDO does not only improve the objective performance, but also reduce its sensitivity to the uncertainty. Therefore, considering the influence of the uncertainty, its typical measurement indexes are the mean and the standard deviation of the objective function. Its mathematical expression can be concisely expressed in(1)minμFor/andσFS.Tgi≤0;i=1,2…,nhj=0j=1,2,…,mwhereμ(F) and σ(F) are the mean and the standard deviation of the objective function respectively, gi is the ith inequality constraint, and hj is the jth equality constraint.


### 2.2. Reliability-Based Design Optimization (RBDO)

The core of RBDO is to deal with the influence of the uncertainty on the constraints, and its mathematical expression is as follows:(2)minFS.TPgi≤0≥r;i=1,2…,nμhj=0j=1,2,…,mwhereP(gi≤0) denotes the probability that the nth constraint is less than or equal to 0, r denotes the reliability that the ith constraint is required to reach, and μ(hj) is the mean of thejth equality constraint.


### 2.3. Reliability-Based Robust Design Optimization

Changing the constraint in RDO to the probability constraint, RBRDO takes the influence of the uncertainty on the objective and the constraint into account simultaneously. Its mathematical expression is as follows:(3)minμFor/andσFS.TPgi≤0≥r;i=1,2…,nμhj=0j=1,2,…,m


## 2.1. Robust Design Optimization (RDO)

Compared with the traditional design optimization, RDO considers the interference of external factors on outputs, whose solutions are not easily perturbed by external factors.As shown in Figure1, point 1 is the optimal point of deterministic optimization (DO) while point 2 is the optimal point of RDO. If the influence of uncertain factors is not considered, the objective function value of point 1 is minimum, which should be the optimal solution of objective function. However, when the design variable x perturbs in the range of ±Δx, the perturbation ΔF1 of objective function at point 1 violates constraints, resulting in failure of the product. The perturbation ΔF2 is obviously less than ΔF1, and all the objective outputs are within the constraint. Therefore, the robust optimal solution is a relatively stable optimal solution.Figure 1
Deterministic and robust optimal solutions.RDO does not only improve the objective performance, but also reduce its sensitivity to the uncertainty. Therefore, considering the influence of the uncertainty, its typical measurement indexes are the mean and the standard deviation of the objective function. Its mathematical expression can be concisely expressed in(1)minμFor/andσFS.Tgi≤0;i=1,2…,nhj=0j=1,2,…,mwhereμ(F) and σ(F) are the mean and the standard deviation of the objective function respectively, gi is the ith inequality constraint, and hj is the jth equality constraint.


## 2.2. Reliability-Based Design Optimization (RBDO)

The core of RBDO is to deal with the influence of the uncertainty on the constraints, and its mathematical expression is as follows:(2)minFS.TPgi≤0≥r;i=1,2…,nμhj=0j=1,2,…,mwhereP(gi≤0) denotes the probability that the nth constraint is less than or equal to 0, r denotes the reliability that the ith constraint is required to reach, and μ(hj) is the mean of thejth equality constraint.


## 2.3. Reliability-Based Robust Design Optimization

Changing the constraint in RDO to the probability constraint, RBRDO takes the influence of the uncertainty on the objective and the constraint into account simultaneously. Its mathematical expression is as follows:(3)minμFor/andσFS.TPgi≤0≥r;i=1,2…,nμhj=0j=1,2,…,m


## 3. Polynomial Chaos Expansions for Sensitivity Analysis

### 3.1. Polynomial Chaos Expansions (PCE)

Compared with MC method, PCE [29] can use the less computational cost to quantify the stochastic property of the output under the influence of multiple uncertain parameters, which has more advantages in ship uncertainty-based design optimization than MC method.Constructing random output variables into a surrogate model with the randomness, for a random variableZ(ε), it can be approximated as (4)Zε=∑i=0∞αiΨiεwhereε=ε1,ε2,…,εn are independent n-dimensional random input variables, Z(ε) is the random output variable, αi are coefficients of polynomial chaos expansions terms, and Ψi(ε) is the multidimensional orthogonal polynomial with respect to ε, which satisfies (5)∫ΩΨiεΨjερεdε=0,i≠jIn order to estimate the uncertainty of outputs, (4) is truncated to a certain order d; therefore, d-order polynomial chaos expansions of Z(ε) is (6)Zε∧≈∑i=0P-1αiΨiεwhere the total number of PCE terms isP=n+d!/n!d!. Since n-dimensional random variables ε are mutually independent, the multidimensional orthogonal polynomial can be decomposed into the product of multiple one-dimensional orthogonal polynomials:(7)Ψiε1,ε2,…,εn=∏t=1nφtltεti=0,1,…,P-1;lt=0,1,…,d;∑t=1nlt≤dwhereφtltεt is thelt-order one-dimensional orthogonal polynomial with respect to εt.In summary, in order to do the multidimensional orthogonal polynomial chaos expansions, coefficients of polynomial chaos expansions termsαi and the one-dimensional orthogonal polynomial basis φtltεt need to be obtained. The solving process of the one-dimensional orthogonal polynomial basis of different probability distributions can be seen in [29]. Ψi(ε) can be obtained from (7) after obtaining the one-dimensional orthogonal polynomial. As for coefficients of polynomial chaos expansions terms αi, given a set of P input sample points ε(0),ε(1),….,ε(P-1) and corresponding outputs Z(ε)(0),Z(ε)(1),…,Z(ε)(P-1), αi can be calculated using (8)Hα=Ψ0ε0Ψ1ε0Ψ0ε1Ψ1ε1⋯ΨP-1ε0⋯ΨP-1ε1⋯⋯Ψ0εP-1Ψ1εP-1⋯⋯⋯ΨP-1εP-1α0α1⋯αP-1=Zε0Zε1⋯ZεP-1=ZwhereH is the polynomial coefficient matrix, α is the coefficient matrix, and Z is the output matrix. Meanwhile, after obtaining coefficients, according to the property of the orthogonal polynomial, the mean and the standard deviation of the output are (9)EZ=α0(10)DZ=∑i=1P-1EΨi2αi2Higher statistical moments can be calculated similarly.


### 3.2. Sobol Decomposition and Indices

First, an input random variable is assumed asx=(x1,x2,…,xn) with joint PDF fx(x). Due to the assumption of independent components of x, fx(x)=∏i=1nfxi(xi), where fxi(xi) is the marginal PDF of xi. For the square integrable function Y(x), it can be decomposed into summands of its main effects and interactions as(11)Yx=Y0+∑i=1nYixi+∑i<jYijxi,xj+⋯+Y12⋯nx1,x2,…,xn=Y0+∑s=1n∑i1<i2<⋯<isnYi1i2⋯isxi1,xi2,…xis=Y0+∑Yuxuwhere1≤i1≤i2≤⋯≤is≤n, u=(i1,i2,…,is). xu represents the subvector of x containing only those components of which the indices belong to u.The Sobol’ decomposition is unique under the condition that the integral of the expansion term for any one of independent variables is 0, as in(12)∫Yuxuftxtdxt=0t∈uEquation (12) leads to the orthogonality property:(13)EYuxuYVxV=0u≠vThe square integral of two sides of (11) is shown in (14)∫Y2xfxxdx-Y02=∑i=1n∫Yi2xifxixidxi+∑i<j∫Yij2xi,xjfxixifxjxjdxidxj+⋯+∫Y12⋯n2x1,x2,…,xnfx1x1fx2x2⋯fxnxndx1dx2⋯dxnThe left side of (14) is the variance of the output and the right side is the partial variance of the output of the corresponding variable. Take first-order partial variance as Vi=∫Yi2(xi)fxi(xi)dxi and second-order variance as Vij=∫Yij2(xi,xj)fxi(xi)fxj(xj)dxidxj. Meanwhile, the uniqueness and orthogonality properties allow decomposition of the variance Vof Y(x) as(15)V=∑i=1nVi+∑i<jVij+⋯+V12⋯n=∑s=1n∑i1<i2<⋯<isnVi1i2⋯is=∑VuwhereVu denotes the partial variance. Then the Sobol’ index SIu is defined as(16)SIu=VuV=EYu2xuVwhere∑SIu=1. For first-order indices, SIi, it describes the influence of each variable xi considered separately as main effects, MSI. Second-order indices, SIij, describe the effects of the interaction between the variables xi and xj. Higher-order indices can be obtained similarly. However, it is obviously impossible to consider all higher-order effects. Therefore, in practical applications, all higher effects are often measured by the total effect sensitivity indices, TSI, accounting for its main effect and all interactions with other parameters.


### 3.3. Sobol’ Indices from Polynomial Chaos Expansions

As mentioned in Section3.1, lt=0,1,…d and ∑t=1nlt≤d. Let us define ϑu as a set of l=(l1,l2,…,ln) and(17)ϑu=l:lt>0∀t=1,2,…,nt∈ult=0∀t=1,2,…,nt∉uOnly whent∈u, l is nonzero. Note that ϑu corresponds to the polynomials depending on u; for example, the polynomial of ϑi is represented by Ψl(εi). Therefore, P-1 terms in (6) may now be gathered according to the parameters they depend on:(18)Zε∧≈∑i=0P-1αiΨiε=α0+∑i=1n∑l∈ϑiαlΨlεi+∑1≤i1<i2≤n∑l∈ϑi1i2αlΨlεi1,εi2+⋯+∑1≤i1<i2<⋯<is≤n∑l∈ϑi1i2⋯isαlΨlεi1,εi2,…,εis+⋯+∑l∈ϑ1,2,…kαlΨlε1,ε2,…εnComparing (18) with (11), the summands in the Sobol’ decomposition of Z(ε)∧ straightforwardly read:(19)Yuxu=∑l∈ϑuαlΨlεIt is now easy to derive sensitivity indices from the above representation. According to (10), (11), (16), and (19), PCE-based Sobol’ indices SIu are defined as(20)SIu=∑l∈ϑuαl2EΨl2DPC


### 3.4. Numerical Examples

#### 3.4.1. Polynomial Model

First, consider a simple polynomial function:(21)Y=12n∏i=1n3Xi2+1where the input variableXi(i=1,2,…n) is uniformly distributed on [0,1]. The analytical solution of Sobol’ indices [14] are(22)SIi1i2⋯is=5-s6/5n-1i=1,2,…,nThis numerical example is carried out forn=3; that is, the order of this polynomial model is 6. Therefore, exact values of sensitivity indices are(23)SI1=SI2=SI3=2591=0.2747SI12=SI23=SI13=591=0.0549SI123=191=0.0110In practice, the true form of the model output is unknown, thus, the order of PCE has to be chosen a priori. Thus 4 cases are considered; namely, the expansion orderd is 3,4,5,6, and sensitivity indices under different orders are listed in Table 1.Table 1
Sensitivity indices of the polynomial model under different expansion orders.
Sensitivity indices Analytical solution Expansion order d = 3 d = 4 d = 5 d = 6 value error value error value error value error S I 1 0.2747 0.2752 0.2% 0.2741 0.2% 0.2747 0 0.2747 0 S I 2 0.2747 0.2752 0.2% 0.2741 0.2% 0.2747 0 0.2747 0 S I 3 0.2747 0.2752 0.2% 0.2741 0.2% 0.2747 0 0.2747 0 S I 12 0.0549 0.0556 1.2% 0.0552 0.5% 0.0550 0.1% 0.0549 0 S I 23 0.0549 0.0556 1.2% 0.0552 0.5% 0.0550 0.1% 0.0549 0 S I 13 0.0549 0.0556 1.2% 0.0552 0.5% 0.0550 0.1% 0.0549 0 S I 123 0.0110 0.0091 17.6% 0.0109 1.0% 0.0110 0.1% 0.0110 0As shown in Table1, compared with the analytical solution, using a 3rd order PCE, the solution has already provided a good estimation of sensitivity indices, since the relative error is less than 0.2% for the first-order indices and 1.2% for second-order indices. For third-order indices, the relative error is larger which reaches 17.6%. When expansion order is 4, these relative errors all drop to 1% or less. The relative error obtained with d=5 and the relative order of the first-order indices were infinitely close to 0 where that of second-order and third-order indices decreases down to 0.1%. When the order of PCE reaches 6, the true order of the function and the relative error are 0 and the exact value is obtained.


#### 3.4.2. Sobol’ Function

Then Sobol’ function is considered:(24)Y=∏i=1q4Xi-2+ai1+aiwhere the input variableXi(i=1,2,…,q) is uniformly distributed over [0,1]. The analytical solution of Sobol’ indices [30] is(25)D=∏i=1qDi+1-1,Di=131+ai2,SIi1i2⋯is=1D∏i=1SDiIn this case,a=[1,2,5,10,20,50,100,500] is selected with q=8. At this time, the order of Sobol’ function is 8. If the 8th PCE with 8 variables is performed, there are 12870 coefficients to be calculated, leading to a high calculation cost. In order to reduce the computational burden, a 2nd order PCE with 45 coefficients is used. MSI and TSI are reported in Table 2 and Figure 2.Table 2
Sensitivity indices of Sobol’ function (d=2).
Sensitivity indices Analytical solution PCE(d=2) M S I 1 0.6037 0.6263 M S I 2 0.2683 0.2814 M S I 3 0.0671 0.0711 M S I 4 0.0200 0.0213 M S I 5 0.0055 0.0062 M S I 6 0.0009 0.0012 M S I 7 0.0002 0.0006 M S I 8 0.0000 0.0002  T S I 1 0.6342 0.6538 T S I 2 0.2945 0.3128 T S I 3 0.0756 0.0792 T S I 4 0.0227 0.0231 T S I 5 0.0062 0.0068 T S I 6 0.0011 0.0015 T S I 7 0.0003 0.0007 T S I 8 0.0000 0.0002Figure 2
Sensitivity indices of Sobol’ function (d=2).It can be seen that first 4 variables are the main influence on the output whether fromMSI or TSI. Therefore, the last 4 variables are fixed at their mean value (X5=X6=X7=X8=0.5) while the first 4 variables are still uniformly distributed over [0,1]. This simplified model is used to replace the original function for analysis with different orders and the results are shown in Table 3.Table 3
Sensitivity indices of simplified model of Sobol’ function with 4 variables under different expansion orders.
Sensitivity indices Analytical solution Expansion order d = 3 d = 5 d = 7 value error value error value error S I 1 0.6037 0.6334 4.9% 0.6144 1.8% 0.6085 0.8% S I 2 0.2683 0.2778 3.5% 0.2680 0.1% 0.2683 0 S I 3 0.0671 0.0685 2.1% 0.0657 2.1% 0.0666 0.7% S I 4 0.0200 0.0203 1.5% 0.0194 3.0% 0.0198 1.0%  T S I 1 0.6342 0.6334 1.3% 0.6434 1.5% 0.6412 1.0% T S I 2 0.2945 0.2778 5.7% 0.2931 0.5% 0.2967 0.7% T S I 3 0.0756 0.0685 9.4% 0.0739 2.2% 0.0756 0 T S I 4 0.0227 0.0203 10.6% 0.0221 2.6% 0.0225 0.9%  Unknown coefficientsP 35 126 330Compared with analytical solutions, some conclusions can be drawn.(1) Whend=3, the relative error of sensitivity indices is large and maximum error reaches nearly 11%. However, as soon as the order of PCE is greater than 5, sensitivity indices have gradually converged to the exact values. Specifically, relative errors of MSI and TSI are both less than 3% for d=5 and 1% for d=7, which shows that the simplified model is able to replace the original model with guaranteed accuracy.(2) For the latter 2 cases, the model evaluations are 171 (45+126) and 375 (45+330), respectively, which is a significant reduction compared with the model of the 8th PCE with 8 variables (12870). Therefore, the simplified model can greatly reduce the computational cost of SA under the premise of ensuring the accuracy; therefore, this method can be applied into practical applications.In conclusion, PCE-based Sobol’ indices have proven the efficiency and accuracy in SA. It is shown that the computation of Sobol’ indices after a proper expansion is precise enough while the computation cost is thus transferred to the obtention of the PCE coefficients, and the subsequent postprocessing being almost costless.


## 3.1. Polynomial Chaos Expansions (PCE)

Compared with MC method, PCE [29] can use the less computational cost to quantify the stochastic property of the output under the influence of multiple uncertain parameters, which has more advantages in ship uncertainty-based design optimization than MC method.Constructing random output variables into a surrogate model with the randomness, for a random variableZ(ε), it can be approximated as (4)Zε=∑i=0∞αiΨiεwhereε=ε1,ε2,…,εn are independent n-dimensional random input variables, Z(ε) is the random output variable, αi are coefficients of polynomial chaos expansions terms, and Ψi(ε) is the multidimensional orthogonal polynomial with respect to ε, which satisfies (5)∫ΩΨiεΨjερεdε=0,i≠jIn order to estimate the uncertainty of outputs, (4) is truncated to a certain order d; therefore, d-order polynomial chaos expansions of Z(ε) is (6)Zε∧≈∑i=0P-1αiΨiεwhere the total number of PCE terms isP=n+d!/n!d!. Since n-dimensional random variables ε are mutually independent, the multidimensional orthogonal polynomial can be decomposed into the product of multiple one-dimensional orthogonal polynomials:(7)Ψiε1,ε2,…,εn=∏t=1nφtltεti=0,1,…,P-1;lt=0,1,…,d;∑t=1nlt≤dwhereφtltεt is thelt-order one-dimensional orthogonal polynomial with respect to εt.In summary, in order to do the multidimensional orthogonal polynomial chaos expansions, coefficients of polynomial chaos expansions termsαi and the one-dimensional orthogonal polynomial basis φtltεt need to be obtained. The solving process of the one-dimensional orthogonal polynomial basis of different probability distributions can be seen in [29]. Ψi(ε) can be obtained from (7) after obtaining the one-dimensional orthogonal polynomial. As for coefficients of polynomial chaos expansions terms αi, given a set of P input sample points ε(0),ε(1),….,ε(P-1) and corresponding outputs Z(ε)(0),Z(ε)(1),…,Z(ε)(P-1), αi can be calculated using (8)Hα=Ψ0ε0Ψ1ε0Ψ0ε1Ψ1ε1⋯ΨP-1ε0⋯ΨP-1ε1⋯⋯Ψ0εP-1Ψ1εP-1⋯⋯⋯ΨP-1εP-1α0α1⋯αP-1=Zε0Zε1⋯ZεP-1=ZwhereH is the polynomial coefficient matrix, α is the coefficient matrix, and Z is the output matrix. Meanwhile, after obtaining coefficients, according to the property of the orthogonal polynomial, the mean and the standard deviation of the output are (9)EZ=α0(10)DZ=∑i=1P-1EΨi2αi2Higher statistical moments can be calculated similarly.


## 3.2. Sobol Decomposition and Indices

First, an input random variable is assumed asx=(x1,x2,…,xn) with joint PDF fx(x). Due to the assumption of independent components of x, fx(x)=∏i=1nfxi(xi), where fxi(xi) is the marginal PDF of xi. For the square integrable function Y(x), it can be decomposed into summands of its main effects and interactions as(11)Yx=Y0+∑i=1nYixi+∑i<jYijxi,xj+⋯+Y12⋯nx1,x2,…,xn=Y0+∑s=1n∑i1<i2<⋯<isnYi1i2⋯isxi1,xi2,…xis=Y0+∑Yuxuwhere1≤i1≤i2≤⋯≤is≤n, u=(i1,i2,…,is). xu represents the subvector of x containing only those components of which the indices belong to u.The Sobol’ decomposition is unique under the condition that the integral of the expansion term for any one of independent variables is 0, as in(12)∫Yuxuftxtdxt=0t∈uEquation (12) leads to the orthogonality property:(13)EYuxuYVxV=0u≠vThe square integral of two sides of (11) is shown in (14)∫Y2xfxxdx-Y02=∑i=1n∫Yi2xifxixidxi+∑i<j∫Yij2xi,xjfxixifxjxjdxidxj+⋯+∫Y12⋯n2x1,x2,…,xnfx1x1fx2x2⋯fxnxndx1dx2⋯dxnThe left side of (14) is the variance of the output and the right side is the partial variance of the output of the corresponding variable. Take first-order partial variance as Vi=∫Yi2(xi)fxi(xi)dxi and second-order variance as Vij=∫Yij2(xi,xj)fxi(xi)fxj(xj)dxidxj. Meanwhile, the uniqueness and orthogonality properties allow decomposition of the variance Vof Y(x) as(15)V=∑i=1nVi+∑i<jVij+⋯+V12⋯n=∑s=1n∑i1<i2<⋯<isnVi1i2⋯is=∑VuwhereVu denotes the partial variance. Then the Sobol’ index SIu is defined as(16)SIu=VuV=EYu2xuVwhere∑SIu=1. For first-order indices, SIi, it describes the influence of each variable xi considered separately as main effects, MSI. Second-order indices, SIij, describe the effects of the interaction between the variables xi and xj. Higher-order indices can be obtained similarly. However, it is obviously impossible to consider all higher-order effects. Therefore, in practical applications, all higher effects are often measured by the total effect sensitivity indices, TSI, accounting for its main effect and all interactions with other parameters.


## 3.3. Sobol’ Indices from Polynomial Chaos Expansions

As mentioned in Section3.1, lt=0,1,…d and ∑t=1nlt≤d. Let us define ϑu as a set of l=(l1,l2,…,ln) and(17)ϑu=l:lt>0∀t=1,2,…,nt∈ult=0∀t=1,2,…,nt∉uOnly whent∈u, l is nonzero. Note that ϑu corresponds to the polynomials depending on u; for example, the polynomial of ϑi is represented by Ψl(εi). Therefore, P-1 terms in (6) may now be gathered according to the parameters they depend on:(18)Zε∧≈∑i=0P-1αiΨiε=α0+∑i=1n∑l∈ϑiαlΨlεi+∑1≤i1<i2≤n∑l∈ϑi1i2αlΨlεi1,εi2+⋯+∑1≤i1<i2<⋯<is≤n∑l∈ϑi1i2⋯isαlΨlεi1,εi2,…,εis+⋯+∑l∈ϑ1,2,…kαlΨlε1,ε2,…εnComparing (18) with (11), the summands in the Sobol’ decomposition of Z(ε)∧ straightforwardly read:(19)Yuxu=∑l∈ϑuαlΨlεIt is now easy to derive sensitivity indices from the above representation. According to (10), (11), (16), and (19), PCE-based Sobol’ indices SIu are defined as(20)SIu=∑l∈ϑuαl2EΨl2DPC


## 3.4. Numerical Examples

### 3.4.1. Polynomial Model

First, consider a simple polynomial function:(21)Y=12n∏i=1n3Xi2+1where the input variableXi(i=1,2,…n) is uniformly distributed on [0,1]. The analytical solution of Sobol’ indices [14] are(22)SIi1i2⋯is=5-s6/5n-1i=1,2,…,nThis numerical example is carried out forn=3; that is, the order of this polynomial model is 6. Therefore, exact values of sensitivity indices are(23)SI1=SI2=SI3=2591=0.2747SI12=SI23=SI13=591=0.0549SI123=191=0.0110In practice, the true form of the model output is unknown, thus, the order of PCE has to be chosen a priori. Thus 4 cases are considered; namely, the expansion orderd is 3,4,5,6, and sensitivity indices under different orders are listed in Table 1.Table 1
Sensitivity indices of the polynomial model under different expansion orders.
Sensitivity indices Analytical solution Expansion order d = 3 d = 4 d = 5 d = 6 value error value error value error value error S I 1 0.2747 0.2752 0.2% 0.2741 0.2% 0.2747 0 0.2747 0 S I 2 0.2747 0.2752 0.2% 0.2741 0.2% 0.2747 0 0.2747 0 S I 3 0.2747 0.2752 0.2% 0.2741 0.2% 0.2747 0 0.2747 0 S I 12 0.0549 0.0556 1.2% 0.0552 0.5% 0.0550 0.1% 0.0549 0 S I 23 0.0549 0.0556 1.2% 0.0552 0.5% 0.0550 0.1% 0.0549 0 S I 13 0.0549 0.0556 1.2% 0.0552 0.5% 0.0550 0.1% 0.0549 0 S I 123 0.0110 0.0091 17.6% 0.0109 1.0% 0.0110 0.1% 0.0110 0As shown in Table1, compared with the analytical solution, using a 3rd order PCE, the solution has already provided a good estimation of sensitivity indices, since the relative error is less than 0.2% for the first-order indices and 1.2% for second-order indices. For third-order indices, the relative error is larger which reaches 17.6%. When expansion order is 4, these relative errors all drop to 1% or less. The relative error obtained with d=5 and the relative order of the first-order indices were infinitely close to 0 where that of second-order and third-order indices decreases down to 0.1%. When the order of PCE reaches 6, the true order of the function and the relative error are 0 and the exact value is obtained.


### 3.4.2. Sobol’ Function

Then Sobol’ function is considered:(24)Y=∏i=1q4Xi-2+ai1+aiwhere the input variableXi(i=1,2,…,q) is uniformly distributed over [0,1]. The analytical solution of Sobol’ indices [30] is(25)D=∏i=1qDi+1-1,Di=131+ai2,SIi1i2⋯is=1D∏i=1SDiIn this case,a=[1,2,5,10,20,50,100,500] is selected with q=8. At this time, the order of Sobol’ function is 8. If the 8th PCE with 8 variables is performed, there are 12870 coefficients to be calculated, leading to a high calculation cost. In order to reduce the computational burden, a 2nd order PCE with 45 coefficients is used. MSI and TSI are reported in Table 2 and Figure 2.Table 2
Sensitivity indices of Sobol’ function (d=2).
Sensitivity indices Analytical solution PCE(d=2) M S I 1 0.6037 0.6263 M S I 2 0.2683 0.2814 M S I 3 0.0671 0.0711 M S I 4 0.0200 0.0213 M S I 5 0.0055 0.0062 M S I 6 0.0009 0.0012 M S I 7 0.0002 0.0006 M S I 8 0.0000 0.0002  T S I 1 0.6342 0.6538 T S I 2 0.2945 0.3128 T S I 3 0.0756 0.0792 T S I 4 0.0227 0.0231 T S I 5 0.0062 0.0068 T S I 6 0.0011 0.0015 T S I 7 0.0003 0.0007 T S I 8 0.0000 0.0002Figure 2
Sensitivity indices of Sobol’ function (d=2).It can be seen that first 4 variables are the main influence on the output whether fromMSI or TSI. Therefore, the last 4 variables are fixed at their mean value (X5=X6=X7=X8=0.5) while the first 4 variables are still uniformly distributed over [0,1]. This simplified model is used to replace the original function for analysis with different orders and the results are shown in Table 3.Table 3
Sensitivity indices of simplified model of Sobol’ function with 4 variables under different expansion orders.
Sensitivity indices Analytical solution Expansion order d = 3 d = 5 d = 7 value error value error value error S I 1 0.6037 0.6334 4.9% 0.6144 1.8% 0.6085 0.8% S I 2 0.2683 0.2778 3.5% 0.2680 0.1% 0.2683 0 S I 3 0.0671 0.0685 2.1% 0.0657 2.1% 0.0666 0.7% S I 4 0.0200 0.0203 1.5% 0.0194 3.0% 0.0198 1.0%  T S I 1 0.6342 0.6334 1.3% 0.6434 1.5% 0.6412 1.0% T S I 2 0.2945 0.2778 5.7% 0.2931 0.5% 0.2967 0.7% T S I 3 0.0756 0.0685 9.4% 0.0739 2.2% 0.0756 0 T S I 4 0.0227 0.0203 10.6% 0.0221 2.6% 0.0225 0.9%  Unknown coefficientsP 35 126 330Compared with analytical solutions, some conclusions can be drawn.(1) Whend=3, the relative error of sensitivity indices is large and maximum error reaches nearly 11%. However, as soon as the order of PCE is greater than 5, sensitivity indices have gradually converged to the exact values. Specifically, relative errors of MSI and TSI are both less than 3% for d=5 and 1% for d=7, which shows that the simplified model is able to replace the original model with guaranteed accuracy.(2) For the latter 2 cases, the model evaluations are 171 (45+126) and 375 (45+330), respectively, which is a significant reduction compared with the model of the 8th PCE with 8 variables (12870). Therefore, the simplified model can greatly reduce the computational cost of SA under the premise of ensuring the accuracy; therefore, this method can be applied into practical applications.In conclusion, PCE-based Sobol’ indices have proven the efficiency and accuracy in SA. It is shown that the computation of Sobol’ indices after a proper expansion is precise enough while the computation cost is thus transferred to the obtention of the PCE coefficients, and the subsequent postprocessing being almost costless.


## 3.4.1. Polynomial Model

First, consider a simple polynomial function:(21)Y=12n∏i=1n3Xi2+1where the input variableXi(i=1,2,…n) is uniformly distributed on [0,1]. The analytical solution of Sobol’ indices [14] are(22)SIi1i2⋯is=5-s6/5n-1i=1,2,…,nThis numerical example is carried out forn=3; that is, the order of this polynomial model is 6. Therefore, exact values of sensitivity indices are(23)SI1=SI2=SI3=2591=0.2747SI12=SI23=SI13=591=0.0549SI123=191=0.0110In practice, the true form of the model output is unknown, thus, the order of PCE has to be chosen a priori. Thus 4 cases are considered; namely, the expansion orderd is 3,4,5,6, and sensitivity indices under different orders are listed in Table 1.Table 1
Sensitivity indices of the polynomial model under different expansion orders.
Sensitivity indices Analytical solution Expansion order d = 3 d = 4 d = 5 d = 6 value error value error value error value error S I 1 0.2747 0.2752 0.2% 0.2741 0.2% 0.2747 0 0.2747 0 S I 2 0.2747 0.2752 0.2% 0.2741 0.2% 0.2747 0 0.2747 0 S I 3 0.2747 0.2752 0.2% 0.2741 0.2% 0.2747 0 0.2747 0 S I 12 0.0549 0.0556 1.2% 0.0552 0.5% 0.0550 0.1% 0.0549 0 S I 23 0.0549 0.0556 1.2% 0.0552 0.5% 0.0550 0.1% 0.0549 0 S I 13 0.0549 0.0556 1.2% 0.0552 0.5% 0.0550 0.1% 0.0549 0 S I 123 0.0110 0.0091 17.6% 0.0109 1.0% 0.0110 0.1% 0.0110 0As shown in Table1, compared with the analytical solution, using a 3rd order PCE, the solution has already provided a good estimation of sensitivity indices, since the relative error is less than 0.2% for the first-order indices and 1.2% for second-order indices. For third-order indices, the relative error is larger which reaches 17.6%. When expansion order is 4, these relative errors all drop to 1% or less. The relative error obtained with d=5 and the relative order of the first-order indices were infinitely close to 0 where that of second-order and third-order indices decreases down to 0.1%. When the order of PCE reaches 6, the true order of the function and the relative error are 0 and the exact value is obtained.


## 3.4.2. Sobol’ Function

Then Sobol’ function is considered:(24)Y=∏i=1q4Xi-2+ai1+aiwhere the input variableXi(i=1,2,…,q) is uniformly distributed over [0,1]. The analytical solution of Sobol’ indices [30] is(25)D=∏i=1qDi+1-1,Di=131+ai2,SIi1i2⋯is=1D∏i=1SDiIn this case,a=[1,2,5,10,20,50,100,500] is selected with q=8. At this time, the order of Sobol’ function is 8. If the 8th PCE with 8 variables is performed, there are 12870 coefficients to be calculated, leading to a high calculation cost. In order to reduce the computational burden, a 2nd order PCE with 45 coefficients is used. MSI and TSI are reported in Table 2 and Figure 2.Table 2
Sensitivity indices of Sobol’ function (d=2).
Sensitivity indices Analytical solution PCE(d=2) M S I 1 0.6037 0.6263 M S I 2 0.2683 0.2814 M S I 3 0.0671 0.0711 M S I 4 0.0200 0.0213 M S I 5 0.0055 0.0062 M S I 6 0.0009 0.0012 M S I 7 0.0002 0.0006 M S I 8 0.0000 0.0002  T S I 1 0.6342 0.6538 T S I 2 0.2945 0.3128 T S I 3 0.0756 0.0792 T S I 4 0.0227 0.0231 T S I 5 0.0062 0.0068 T S I 6 0.0011 0.0015 T S I 7 0.0003 0.0007 T S I 8 0.0000 0.0002Figure 2
Sensitivity indices of Sobol’ function (d=2).It can be seen that first 4 variables are the main influence on the output whether fromMSI or TSI. Therefore, the last 4 variables are fixed at their mean value (X5=X6=X7=X8=0.5) while the first 4 variables are still uniformly distributed over [0,1]. This simplified model is used to replace the original function for analysis with different orders and the results are shown in Table 3.Table 3
Sensitivity indices of simplified model of Sobol’ function with 4 variables under different expansion orders.
Sensitivity indices Analytical solution Expansion order d = 3 d = 5 d = 7 value error value error value error S I 1 0.6037 0.6334 4.9% 0.6144 1.8% 0.6085 0.8% S I 2 0.2683 0.2778 3.5% 0.2680 0.1% 0.2683 0 S I 3 0.0671 0.0685 2.1% 0.0657 2.1% 0.0666 0.7% S I 4 0.0200 0.0203 1.5% 0.0194 3.0% 0.0198 1.0%  T S I 1 0.6342 0.6334 1.3% 0.6434 1.5% 0.6412 1.0% T S I 2 0.2945 0.2778 5.7% 0.2931 0.5% 0.2967 0.7% T S I 3 0.0756 0.0685 9.4% 0.0739 2.2% 0.0756 0 T S I 4 0.0227 0.0203 10.6% 0.0221 2.6% 0.0225 0.9%  Unknown coefficientsP 35 126 330Compared with analytical solutions, some conclusions can be drawn.(1) Whend=3, the relative error of sensitivity indices is large and maximum error reaches nearly 11%. However, as soon as the order of PCE is greater than 5, sensitivity indices have gradually converged to the exact values. Specifically, relative errors of MSI and TSI are both less than 3% for d=5 and 1% for d=7, which shows that the simplified model is able to replace the original model with guaranteed accuracy.(2) For the latter 2 cases, the model evaluations are 171 (45+126) and 375 (45+330), respectively, which is a significant reduction compared with the model of the 8th PCE with 8 variables (12870). Therefore, the simplified model can greatly reduce the computational cost of SA under the premise of ensuring the accuracy; therefore, this method can be applied into practical applications.In conclusion, PCE-based Sobol’ indices have proven the efficiency and accuracy in SA. It is shown that the computation of Sobol’ indices after a proper expansion is precise enough while the computation cost is thus transferred to the obtention of the PCE coefficients, and the subsequent postprocessing being almost costless.


## 4. Probabilistic Collocation Method (PCM)

As mentioned in Section3, a set of sampling points need to be generated to solve the coefficient of PCE. The number of sampling points determines the calculation times of outputs, which directly affects the calculation cost. In addition, the distribution of sampling points will also affect the fitting effect of PCE; thus this affects the accuracy of subsequent results (statistical moments). Due to the orthogonality of PCE, probabilistic collocation method (PCM) is adopted instead of commonly used statistical sampling methods (LHS, etc.)


### 4.1. Implementation of PCM

Consideringε as one-dimensional variable, the residual between the truncated model Z(ε)∧ and the actual model Z(ε) is defined as(26)Rα,ε=Zε-Zε∧The coefficientαi can be calculated by the weighted residual method; therefore, the residual R(αi,ε) and the one-dimensional orthogonal polynomial φi(ε) should be orthogonal to each other:(27)∫ΩRα,εφiεdε=0Equation (12) can be solved by using the Gaussian quadrature approximation:(28)∫ΩRα,εφiεdε≈∑j=0P-1wjRα,εjφiεj=0i=0,1,…,P-1;wherewj and εj are weights and abscissas of Gaussian integral points, respectively. If wjφi(εj) has the same sign and is not zero for all i and j, (13) can be simplified as (29)Rα,εj=0j=0,1,…,P-1;Equation (29) is an equation with respect to only coefficients and abscissas of Gaussian integral point, which represents that the coefficient can be solved by using the known Gaussian integral point.Therefore, ford-order polynomial expansions of one-dimensional random variables, collocation points are simply the roots of one-dimensional orthogonal polynomial with the order of d+1. For instance, Hermite polynomials are the optimal polynomial of the normal distribution, which has the highest probability at the origin. Therefore, collocation points should be arranged as close and symmetrical as possible to the origin. In addition, odd-order PCE requires roots of even-order Hermite polynomials as collocation points, whereas even-order polynomials have no zero root; therefore, zero root is added as a collocation point of odd-order PCE [31]. Thus, the root of the first 4 order Hermite polynomials are shown in Table 4.Table 4
Roots of Hermite polynomials.
Expansion orderd d + 1 order Hermite polynomial Roots (collocation points) 1 ε 2 - 1 0,1 , - 1 2 ε 3 - 3 ε 0 , 3 , - 3 3 ε 4 - 6 ε 2 + 4 0 , 3 + 6 , 3 - 6 , - 3 + 6 , - 3 - 6 4 ε 5 - 10 ε 3 + 15 ε 0 , 5 + 10 , 5 - 10 , - 5 + 10 , - 5 - 10For n-dimensional random variables, collocation points are the combination of roots of the next higher-order orthogonal polynomial. Then the number of all the available probability distribution points is(30)N=d+1n+1,disodd0,diseven


### 4.2. Improvement in PCM

Input samples of PCM are selected according to some certain rules. Generally speaking, the number of sample points should be greater than the number of undetermined coefficients. Isukappali [25, 26] thought that the number of collocation points should be twice as many as the number of undetermined coefficients and then used the regression method to calculate coefficients. Unfortunately, Jiang et al. [27] found that sometimes increasing the number of collocation points to more than 8 times the number of undetermined coefficients still cannot maintain the calculation accuracy. The reason is that the polynomial coefficient matrix contains many linearly dependent collocation points and these points are not helpful to improve the calculation accuracy. Rejecting these linearly dependent points will improve the computational efficiency; therefore, based on the principle of linearly independence, an improved probabilistic collocation method (IPCM) is proposed to reject these linearly dependent points and ensure the linearly independence of the coefficient matrix; that is, the matrix is the full row rank matrix. Figure 3 is the flow chart of this method and the steps are as follows.Figure 3
Flow chart of IPCM.Step 1.
The expansion order, the distribution type, and the number of variables are determined.Step 2.
The collocation points are arranged in the decreasing probability density (an ascending order according their distance from the origin). For instance, the point (0,0,…,0) should be selected first.Step 3.
Then the sorted collocation points are selected one by one to establish the polynomial coefficient matrix row by row.Step 4.
For the candidate of the (i+1)th collocation point, the (i+1)th row of the matrix should be linearly independent with the previous ith rows. Therefore, the rank of this matrix is calculated and judged whether it is equal to the row of this matrix. If it is equal, this candidate is preserved or otherwise abandoned. Then the point with next highest probability density should be tested.Step 5.
The process does not end until the number of the row is equal to the number of undetermined coefficients.When the number of random variables and the degree of PCE are both high, due to the multiple times of composing the polynomial coefficient matrix and calculating the rank of the matrix, the calculation cost of searching linearly independent points can be high as well after following the above steps. Fortunately, if the same basis polynomial chaos is adopted, for a given number of random variables and the order of PCE, collocation points for computing different outputs are the same, which means that linearly independent collocation points can be searched and saved in advance and can be used whenever you need to save the time of repeating searching the same collocation points.


### 4.3. Numerical Examples

To demonstrate the uncertainty quantification capability and illustrate the application of IPCM, a nonlinear nonmonotonicIshigami function is considered, using IPCM and PCM to select collocation points, respectively. (31)Y=sinX1+asin2X2+bX34sinX1where the input variableXi(i=1,2,3) is uniformly distributed on [-π,π] and in this case a=7, b=0.1. Due to its nonlinearity, its true order cannot be determined. Thus 4 cases are considered; namely, the expansion order d is 3,5,7,9 and the number of their corresponding expansion terms P is 20, 56, 120, and 220, respectively. Taking the result of MC method (1000 sample points) as the exact value, PCM (the number of collocation points is 125, 343, 729, and 1331, respectively) and IPCM are used to select sample points for X1 and X2 to calculate the mean and the standard deviation of Y. The variation of them with the expansion order is presented in Figure 4.Figure 4
The variation of the mean and the standard deviation of Y with the expansion order.From Figure4, it can be seen that errors of the mean and the standard deviation calculated by PCM and IPCM are large when the expansion order is low (d=3,5). When the order rises, the mean and the standard deviation tend to converge to the result of MC method. Specifically, when d=7, the mean and the standard deviation calculated by these 2 methods are within 1% error band of the result of MC method. When d=9, the accuracy of the mean and the standard deviation slightly increases compared with d=7; however, the number of collocation points is more than 2 times that of d=7. Therefore, PCE have different optimal expansion orders for different problems; if the expansion order increases blindly, it may obtain a more accurate result but result in higher computational cost.Meanwhile, under the same expansion order, although the accuracy of PCM is slightly higher than that of IPCM, the number of collocation points required by IPCM is much less than that of PCM. As the order increases, for example,d=9, the number of collocation points required by PCM (1331) is even 6 times that of IPCM (220). In conclusion, IPCM has an obvious advantage in the calculation; that is, it can guarantee the calculation accuracy with much fewer collocation points, compared with PCM and MC method.As can be seen, the number of collocation points selected by PCM is greater than the number of undetermined coefficients. However, IPCM can maintain the accuracy by selecting a part of the points from collocation points generated by PCM. In order to further explain the relationship between the number of collocation points and the calculation accuracy, a part of collocation points generated by PCM, namely, 2 times, 3 times, and 4 times the number of undetermined coefficients, is selected to calculate the mean and the standard deviation by the regression method.Tables5 and 6 give the comparison of the mean and the standard deviation of Y for the different number of collocation points. When d=3,5,9, the increase in the number of collocation points does not improve the calculation accuracy; when d=7, the calculation accuracy can only be improved when the number of collocation points increases to 3P.Table 5
Comparison of the mean with the different number of collocation points.
Number of collocation points (N) PCE MC d = 3 d = 5 d = 7 d = 9 (N=1000) N=2P 2.591 (25.6%) 3.098 (11.1%) 2.925 (16.1%) 3.503 (0.5%) N=3P 2.614 (24.9%) 3.112 (10.7%) 3.441 (1.3%) 3.478 (0.2%) N=4P 2.685 (22.9%) 3.159 (9.4%) 3.468 (1.3%) 3.491 (0.2%) 3.485  Full rank 2.912 (19.3%) 3.279 (5.9%) 3.472 (0.4%) 3.489 (0.1%) IPCM 2.676 (23.2%) 3.185 (8.6%) 3.450 (1.0%) 3.495 (0.2%)Table 6
Comparison of standard deviation with different number of collocation points.
Number of collocation points (N) PCE MC d = 3 d = 5 d = 7 d = 9 (N=1000) N=2P 2.786 (24.5%) 3.352 (9.1%) 3.016 (18.3%) 3.743 (1.4%) N=3P 2.828 (23.4%) 3.401 (7.8%) 3.639 (1.4%) 3.735 (1.1%) N=4P 2.953 (19.9%) 3.414 (7.5%) 3.645 (1.2%) 3.728 (1.0%) 3.691  Full rank 3.014 (18.3%) 3.536 (4.2%) 3.727 (0.9%) 3.699 (0.2%) IPCM 2.711 (26.5%) 3.421 (7.3%) 3.728 (1.0%) 3.725 (0.9%)In order to explain the above relationship between the calculation accuracy and the number of collocation points, Figure5 gives the relationship between the rank of the coefficient matrix and the number of collocation points. It can be observed that when the number of collocation points is greater than the number of undetermined coefficients, the polynomial coefficient matrix may not reach the full rank. In this case, errors of results in Tables 5 and 6 are large. As shown in Figure 5, when d=3, the polynomial coefficient matrix can reach the full rank when the number of collocation points is greater than 30. Then the rank of the coefficient matrix remains unchanged even if the number of collocation points increases, which means that increasing the number of collocation points will not improve the calculation accuracy anymore. When N=4P=80, relative errors of the mean and the standard deviation are 22.9% and 19.9%, respectively, while relative errors when N=30 are 19.3% and 18.3%. In these 2 cases, the polynomial coefficient matrices both have already reached the full rank; thus the accuracy does not improve significantly. Nevertheless, the number of collocation points of the former is more than 2 times that of the latter. Similar conclusions can be drawn when d=5,9.Figure 5
Relationship between the rank of the polynomial coefficient matrix and the number of collocation points.However, whend=7, the polynomial coefficient matrix reaches the full rank when the number of collocation points are greater than 247, which explains the improvement in calculation accuracy when N=3P=360. When N=2P, the number of collocation points is less than 247 and the polynomial coefficient matrix does not reach the full rank resulting in large errors even no solution. Besides, the rank of the coefficient matrix does not increase monotonically with the increase of the number of collocation points. Taking the case d=7 as an example, from the 172th to the 235th collocation point, the rank of the polynomial coefficient matrix is always 108, which indicates that these collocation points are linearly dependent and only one point is selected while other points can be eliminated.When using the regression method, the minimum number of collocation points that guarantee the full rank of the coefficient is 30(d=3), 84(d=5), 247(d=7), and 404(d=9). However, the number of collocation points that guarantees the full rank of the coefficient is 20(d=3), 56(d=5), 120(d=7), and 240(d=9) using IPCM. It can be seen that IPCM can reduce the number of collocation points required for the calculation greatly, thereby leading to a decline in the computational cost. The accuracy of the former is slightly higher than that of the latter; this is because the latter only uses a part of collocation points of the former. Nevertheless, compared with the exact solution, the result of IPCM can fully meet the requirement of the accuracy, and the computational efficiency is improved.In summary, when selecting a part of collocation points generated by PCM to solve undetermined coefficients of the polynomial, the selected collocation points should satisfy the condition that the rank of the coefficient matrix should be equal to the number of the undetermined coefficient. However, for the high-dimensional high-order problems, complex uncertainty analyses, and simulations, a large amount of calculation is required. In these cases, IPCM has great advantages due to the fact that IPCM based on the linear independent principle only need the same number of collocation points as the number of undetermined coefficients, with an obvious decline in required collocation points compared with PCM.


## 4.1. Implementation of PCM

Consideringε as one-dimensional variable, the residual between the truncated model Z(ε)∧ and the actual model Z(ε) is defined as(26)Rα,ε=Zε-Zε∧The coefficientαi can be calculated by the weighted residual method; therefore, the residual R(αi,ε) and the one-dimensional orthogonal polynomial φi(ε) should be orthogonal to each other:(27)∫ΩRα,εφiεdε=0Equation (12) can be solved by using the Gaussian quadrature approximation:(28)∫ΩRα,εφiεdε≈∑j=0P-1wjRα,εjφiεj=0i=0,1,…,P-1;wherewj and εj are weights and abscissas of Gaussian integral points, respectively. If wjφi(εj) has the same sign and is not zero for all i and j, (13) can be simplified as (29)Rα,εj=0j=0,1,…,P-1;Equation (29) is an equation with respect to only coefficients and abscissas of Gaussian integral point, which represents that the coefficient can be solved by using the known Gaussian integral point.Therefore, ford-order polynomial expansions of one-dimensional random variables, collocation points are simply the roots of one-dimensional orthogonal polynomial with the order of d+1. For instance, Hermite polynomials are the optimal polynomial of the normal distribution, which has the highest probability at the origin. Therefore, collocation points should be arranged as close and symmetrical as possible to the origin. In addition, odd-order PCE requires roots of even-order Hermite polynomials as collocation points, whereas even-order polynomials have no zero root; therefore, zero root is added as a collocation point of odd-order PCE [31]. Thus, the root of the first 4 order Hermite polynomials are shown in Table 4.Table 4
Roots of Hermite polynomials.
Expansion orderd d + 1 order Hermite polynomial Roots (collocation points) 1 ε 2 - 1 0,1 , - 1 2 ε 3 - 3 ε 0 , 3 , - 3 3 ε 4 - 6 ε 2 + 4 0 , 3 + 6 , 3 - 6 , - 3 + 6 , - 3 - 6 4 ε 5 - 10 ε 3 + 15 ε 0 , 5 + 10 , 5 - 10 , - 5 + 10 , - 5 - 10For n-dimensional random variables, collocation points are the combination of roots of the next higher-order orthogonal polynomial. Then the number of all the available probability distribution points is(30)N=d+1n+1,disodd0,diseven


## 4.2. Improvement in PCM

Input samples of PCM are selected according to some certain rules. Generally speaking, the number of sample points should be greater than the number of undetermined coefficients. Isukappali [25, 26] thought that the number of collocation points should be twice as many as the number of undetermined coefficients and then used the regression method to calculate coefficients. Unfortunately, Jiang et al. [27] found that sometimes increasing the number of collocation points to more than 8 times the number of undetermined coefficients still cannot maintain the calculation accuracy. The reason is that the polynomial coefficient matrix contains many linearly dependent collocation points and these points are not helpful to improve the calculation accuracy. Rejecting these linearly dependent points will improve the computational efficiency; therefore, based on the principle of linearly independence, an improved probabilistic collocation method (IPCM) is proposed to reject these linearly dependent points and ensure the linearly independence of the coefficient matrix; that is, the matrix is the full row rank matrix. Figure 3 is the flow chart of this method and the steps are as follows.Figure 3
Flow chart of IPCM.Step 1.
The expansion order, the distribution type, and the number of variables are determined.Step 2.
The collocation points are arranged in the decreasing probability density (an ascending order according their distance from the origin). For instance, the point (0,0,…,0) should be selected first.Step 3.
Then the sorted collocation points are selected one by one to establish the polynomial coefficient matrix row by row.Step 4.
For the candidate of the (i+1)th collocation point, the (i+1)th row of the matrix should be linearly independent with the previous ith rows. Therefore, the rank of this matrix is calculated and judged whether it is equal to the row of this matrix. If it is equal, this candidate is preserved or otherwise abandoned. Then the point with next highest probability density should be tested.Step 5.
The process does not end until the number of the row is equal to the number of undetermined coefficients.When the number of random variables and the degree of PCE are both high, due to the multiple times of composing the polynomial coefficient matrix and calculating the rank of the matrix, the calculation cost of searching linearly independent points can be high as well after following the above steps. Fortunately, if the same basis polynomial chaos is adopted, for a given number of random variables and the order of PCE, collocation points for computing different outputs are the same, which means that linearly independent collocation points can be searched and saved in advance and can be used whenever you need to save the time of repeating searching the same collocation points.


## 4.3. Numerical Examples

To demonstrate the uncertainty quantification capability and illustrate the application of IPCM, a nonlinear nonmonotonicIshigami function is considered, using IPCM and PCM to select collocation points, respectively. (31)Y=sinX1+asin2X2+bX34sinX1where the input variableXi(i=1,2,3) is uniformly distributed on [-π,π] and in this case a=7, b=0.1. Due to its nonlinearity, its true order cannot be determined. Thus 4 cases are considered; namely, the expansion order d is 3,5,7,9 and the number of their corresponding expansion terms P is 20, 56, 120, and 220, respectively. Taking the result of MC method (1000 sample points) as the exact value, PCM (the number of collocation points is 125, 343, 729, and 1331, respectively) and IPCM are used to select sample points for X1 and X2 to calculate the mean and the standard deviation of Y. The variation of them with the expansion order is presented in Figure 4.Figure 4
The variation of the mean and the standard deviation of Y with the expansion order.From Figure4, it can be seen that errors of the mean and the standard deviation calculated by PCM and IPCM are large when the expansion order is low (d=3,5). When the order rises, the mean and the standard deviation tend to converge to the result of MC method. Specifically, when d=7, the mean and the standard deviation calculated by these 2 methods are within 1% error band of the result of MC method. When d=9, the accuracy of the mean and the standard deviation slightly increases compared with d=7; however, the number of collocation points is more than 2 times that of d=7. Therefore, PCE have different optimal expansion orders for different problems; if the expansion order increases blindly, it may obtain a more accurate result but result in higher computational cost.Meanwhile, under the same expansion order, although the accuracy of PCM is slightly higher than that of IPCM, the number of collocation points required by IPCM is much less than that of PCM. As the order increases, for example,d=9, the number of collocation points required by PCM (1331) is even 6 times that of IPCM (220). In conclusion, IPCM has an obvious advantage in the calculation; that is, it can guarantee the calculation accuracy with much fewer collocation points, compared with PCM and MC method.As can be seen, the number of collocation points selected by PCM is greater than the number of undetermined coefficients. However, IPCM can maintain the accuracy by selecting a part of the points from collocation points generated by PCM. In order to further explain the relationship between the number of collocation points and the calculation accuracy, a part of collocation points generated by PCM, namely, 2 times, 3 times, and 4 times the number of undetermined coefficients, is selected to calculate the mean and the standard deviation by the regression method.Tables5 and 6 give the comparison of the mean and the standard deviation of Y for the different number of collocation points. When d=3,5,9, the increase in the number of collocation points does not improve the calculation accuracy; when d=7, the calculation accuracy can only be improved when the number of collocation points increases to 3P.Table 5
Comparison of the mean with the different number of collocation points.
Number of collocation points (N) PCE MC d = 3 d = 5 d = 7 d = 9 (N=1000) N=2P 2.591 (25.6%) 3.098 (11.1%) 2.925 (16.1%) 3.503 (0.5%) N=3P 2.614 (24.9%) 3.112 (10.7%) 3.441 (1.3%) 3.478 (0.2%) N=4P 2.685 (22.9%) 3.159 (9.4%) 3.468 (1.3%) 3.491 (0.2%) 3.485  Full rank 2.912 (19.3%) 3.279 (5.9%) 3.472 (0.4%) 3.489 (0.1%) IPCM 2.676 (23.2%) 3.185 (8.6%) 3.450 (1.0%) 3.495 (0.2%)Table 6
Comparison of standard deviation with different number of collocation points.
Number of collocation points (N) PCE MC d = 3 d = 5 d = 7 d = 9 (N=1000) N=2P 2.786 (24.5%) 3.352 (9.1%) 3.016 (18.3%) 3.743 (1.4%) N=3P 2.828 (23.4%) 3.401 (7.8%) 3.639 (1.4%) 3.735 (1.1%) N=4P 2.953 (19.9%) 3.414 (7.5%) 3.645 (1.2%) 3.728 (1.0%) 3.691  Full rank 3.014 (18.3%) 3.536 (4.2%) 3.727 (0.9%) 3.699 (0.2%) IPCM 2.711 (26.5%) 3.421 (7.3%) 3.728 (1.0%) 3.725 (0.9%)In order to explain the above relationship between the calculation accuracy and the number of collocation points, Figure5 gives the relationship between the rank of the coefficient matrix and the number of collocation points. It can be observed that when the number of collocation points is greater than the number of undetermined coefficients, the polynomial coefficient matrix may not reach the full rank. In this case, errors of results in Tables 5 and 6 are large. As shown in Figure 5, when d=3, the polynomial coefficient matrix can reach the full rank when the number of collocation points is greater than 30. Then the rank of the coefficient matrix remains unchanged even if the number of collocation points increases, which means that increasing the number of collocation points will not improve the calculation accuracy anymore. When N=4P=80, relative errors of the mean and the standard deviation are 22.9% and 19.9%, respectively, while relative errors when N=30 are 19.3% and 18.3%. In these 2 cases, the polynomial coefficient matrices both have already reached the full rank; thus the accuracy does not improve significantly. Nevertheless, the number of collocation points of the former is more than 2 times that of the latter. Similar conclusions can be drawn when d=5,9.Figure 5
Relationship between the rank of the polynomial coefficient matrix and the number of collocation points.However, whend=7, the polynomial coefficient matrix reaches the full rank when the number of collocation points are greater than 247, which explains the improvement in calculation accuracy when N=3P=360. When N=2P, the number of collocation points is less than 247 and the polynomial coefficient matrix does not reach the full rank resulting in large errors even no solution. Besides, the rank of the coefficient matrix does not increase monotonically with the increase of the number of collocation points. Taking the case d=7 as an example, from the 172th to the 235th collocation point, the rank of the polynomial coefficient matrix is always 108, which indicates that these collocation points are linearly dependent and only one point is selected while other points can be eliminated.When using the regression method, the minimum number of collocation points that guarantee the full rank of the coefficient is 30(d=3), 84(d=5), 247(d=7), and 404(d=9). However, the number of collocation points that guarantees the full rank of the coefficient is 20(d=3), 56(d=5), 120(d=7), and 240(d=9) using IPCM. It can be seen that IPCM can reduce the number of collocation points required for the calculation greatly, thereby leading to a decline in the computational cost. The accuracy of the former is slightly higher than that of the latter; this is because the latter only uses a part of collocation points of the former. Nevertheless, compared with the exact solution, the result of IPCM can fully meet the requirement of the accuracy, and the computational efficiency is improved.In summary, when selecting a part of collocation points generated by PCM to solve undetermined coefficients of the polynomial, the selected collocation points should satisfy the condition that the rank of the coefficient matrix should be equal to the number of the undetermined coefficient. However, for the high-dimensional high-order problems, complex uncertainty analyses, and simulations, a large amount of calculation is required. In these cases, IPCM has great advantages due to the fact that IPCM based on the linear independent principle only need the same number of collocation points as the number of undetermined coefficients, with an obvious decline in required collocation points compared with PCM.


## 5. Preliminary Design of the Bulk Carrier

### 5.1. Description of Optimization Design

The proposed method is applied to the ship UDO. The optimization model of the bulk carrier proposed by Hannapel S and Vlahopoulos N is considered. Details of this model can be found in [32]. Deterministic optimization (DO) and uncertainty optimization are both applied to this model.


### 5.2. Specific Process of Bulk Carrier Uncertainty-Based Design Optimization

The specific process is shown in Figure6. The steps of this process are as follows.Figure 6
Flow chart of ship uncertainty-based design optimization.Step 1.
Preparation work: determine the optimal PCE order of original model, then perform the sensitivity analysis on uncertain variables and rank them according to their degree of influence on the objective and the constraint, fix the uncertain variable with less influence at the mean value, and reduce the dimension of the model. Redetermine the optimal PCE order for the simplified model.Step 2.
Design variables and their ranges are determined.Step 3.
A set of design variables, in other words, a bulk carrier scheme, is selected and transferred to the inner layer, the uncertainty analysis module.Step 4.
According to the expansion order, the distribution type, and the number of uncertain variables, which are determined in advance, IPCM is used to select collocation points for uncertain variables.Step 5.
PCE is used to perform the expansion of the objective and the constraint.Step 6.
According to the design requirement, select one of 3 types of UDO. Then the mean and the standard deviation of the objective and the failure probability of the constraint are obtained and transferred to the outer layer, optimization module.Step 7.
The optimizer judges whether the entire optimization process meets requirements of the iteration number. If the process reaches the iteration number, the optimization process ends. If not, step 3~step 6 are implemented again.


### 5.3. Optimization Analyses

#### 5.3.1. Deterministic Optimization

In order to give a baseline for the comparison for later results of uncertainty optimization, DO is done first. DO uses a multi-island genetic algorithm, with the subpopulation size being 10, the number of islands being 10, and the number of generations being 50. Design variables, constraints, and objective functions are presented in Table7 and the comparison between results of DO and initial design [32] is listed in Table 12.Table 7
Bulk carrier model.
Design variables L e n g t h , L / m 150 ≤ L ≤ 274.32 B r e a d t h , B / m 20 ≤ B ≤ 32.31 D e p t h , D / m 13 ≤ D ≤ 25 D r a f t , T / m 10 ≤ D ≤ 11.71 B l o c k   c o e f f i c i e n t , C B 0.63 ≤ C B ≤ 0.75 D e s i g n   s p e e d , V k 14 ≤ V k ≤ 20  Objective function  T r a n s p o r t a t i o n   c o s t , T C / ( P o u n d / y r ) T C = C A / A C  constraints  L / B ≥ 6 L / D ≤ 15 L / T ≤ 19 F n ≤ 0.32 25000 ≤ D W ≤ 500000 T - 0.45 D W 0.31 ≤ 0 T - 0.7 D - 0.7 ≤ 0 0.07 B - 0.35 T - ( 0.085 C B - 0.002 ) B 2 T C B + 1 + 0.25 D ≤ 0


#### 5.3.2. Uncertainty-Based Optimization

Selection of Uncertain Variables. In the real travelling process, due to the unavoidable disturbance such as the wind and the wave, the ship resistance in waves will increase with more violent ship motions; therefore the service/operational speed of ship will decrease compared with the hydrostatic condition. Therefore, the speed of a sailing ship is not a fixed value, but an uncertain variable that perturbs around the design speed. By monitoring the ship speed in a period of time [6, 13], we obviously find that the speed perturbation is a universal existence. Although the perturbation uncertainty magnitudes are relatively small in most cases, their continuous calculation and coupling effects with other parameters will cause the system response to a large deviation in the optimization design process. Therefore, the ship speed is selected as an uncertain variable in (32)Vactual=Vk+ZwhereVactual is the actual speed, Vk is the design speed, and Z is the speed perturbation, which follows the normal distribution with the mean of 0 and the standard deviation of 0.1Vk [6].Then the uncertainty in parameters is introduced. The steel weight of the bulk carrier is given by(33)WS=0.034L1.7B0.7D0.4CB0.5Equation (33) is obtained by the statistical regression using the real ship data. Considering the limited number of samples and the uncontrollable error of using the unknown regression method, the regression exponent cannot accurately reflect the actual engineering situation; therefore, the uncertainty is applied to the exponent [32], which are assumed to be random parameters, namely, X (the normal distribution with the mean of 1.7 and the standard deviation of 0.17) and Y (the normal distribution with the mean of 0.7 and the standard deviation of 0.07), A (the normal distribution with the mean of 0.4 and the standard deviation of 0.04), and B (the normal distribution with the mean of 0.5 and the standard deviation of 0.05), as presented in (34)WS=0.034LXBYDACBBDetermination of the Optimal PCE Order. Considering the influence of the uncertainty on the objective and the constraint for the DO optimal case, the objective and the constraint are computed using PCE, whose expansion order is d=2,3,4 with the number of their corresponding expansion terms P are 21, 56, and 126, respectively. IPCM is used to select collocation points, for the optimization objective, and the constraint, the value of undetermined coefficients of order d=2,3,4 is shown in Figure 7, where the abscissa is the serial number of the coefficient and the ordinate is the coefficient value.Figure 7
Coefficient values of the objective and the constraint with different expansion orders.It can be seen from Figure10 that, for T-0.45DW0.31, undetermined coefficient values after the number 72 are almost all 0; thus the 4th order expansion can basically meet the requirement of the accuracy. For Fn, the coefficient values of d=2,3,4 are nearly the same, thus the 2nd-order expansion is optimal for Fn. For DW, the coefficient values after the number 56 (the 4th terms) are almost all 0, which means the 4th terms have no contribution to DW; therefore, 3rd-order expansion is selected. For TC, coefficient values have been close to 0 since the number 107; therefore, the 4th-order can basically meet requirement of the accuracy. In summary, the optimal PCE order is the 4th-order.Sensitivity Analysis. Then, for the objective and the constraint, sensitivity analysis is performed by using PCE method and traditional Sobol’ method. The sensitivity indices, namely, MSI and TSI of the objective and the constraint to X, Y, Z, A, and B, are shown in Table 8 and Figure 8.Table 8
Sensitivity indices of the objective and constraints.
T - 0.45 D W 0.31 D W F n T C PCE Sobol’ PCE Sobol’ PCE Sobol PCE Sobol’ M S I X 0.8692 0.8857 0.8842 0.8302 0 0 0.8507 0.8925 M S I Y 0.0615 0.0637 0.0595 0.0582 0 0 0.0649 0.0712 M S I Z 0.0002 0.0012 0.0002 0.0004 0.9999 0.9249 0.0023 0.0021 M S I A 0.0136 0.0153 0.0128 0.0139 0 0 0.0148 0.151 M S I B 0.0006 0.0005 0.0005 0.0005 0 0 0.0007 0.0006 T S I X 0.9238 0.9159 0.9267 0.9451 0 0 0.9169 0.9237 T S I Y 0.1060 0.1155 0.0945 0.1151 0 0 0.1180 0.1244 T S I Z 0.0002 0.0003 0.0002 0.0002 0.9999 0.9249 0.0027 0.0031 T S I A 0.0241 0.0230 0.0208 0.0314 0 0 0.0279 0.0293 T S I B 0.0011 0.0020 0.0009 0.0012 0 0 0.0013 0.0021Figure 8
Sensitivity indices of the objective and constraints.According toMSI and TSI, 5 uncertain variables are ranked in order. For T-0.45DW0.31, DW and TC, X, and Y have dominated the role and A has only a little effect, while the influence of other variables can be neglected. On the contrary, for Fn, Z appears to be the most important in the decomposition of the variance of the response while the effect of other 4 variables has no need to consider, which is consistent with the formula of Fn. In conclusion, whether from MSI or TSI, X, Y, and Z have main effects on the objective and the constraint. Therefore, A and B are fixed at their mean values, while X, Y, and Z still satisfy their original distributions. The simplified model is used to replace the original model to perform optimization.For this simplified model, the process of determining the optimal PCE order is repeated, and 3rd-order expansion is optimal for the model. The mean and the standard deviation of 3rd-order expansion are compared with results of MC method, as can be seen in Table9.Table 9
The mean and the standard deviation of the objective and the constraint.
the objective /the constraint T C F n ≤ 0.32 25000 ≤ D W ≤ 500000 T - 0.45 D W 0.31 ≤ 0 Mean MC (100) 10.06 0.1699 3.408∗104 0.309 PCE (20) 9.92 0.1699 3.479∗104 0.304 Error 1.39 % 0 % 1.95 % 1.64 %  Standard deviation MC (100) 4.21 0.017 8009 0.988 PCE (20) 4.29 0.017 7968 1.021 Error 1.90 % 0 % 0.52 % 3.18 %From results in Table9, taking the result of MC method as the exact value, the error of the mean and the standard deviation solved by PCE is within 2% and 3.5%, respectively. Meanwhile required sample points of PCE (20 sample points) are far less than those of MC method (100 sample points). Therefore, the advantage of PCE in computational efforts is obvious; specifically, required sample points reduce greatly while the calculation accuracy is not significantly reduced.Then, according to the solved first 4 order moments of the constraints, maximum entropy method (MEM) is adopted to obtain PDF of constraints and then the failure probability is solved. Compared with MC method, the results are listed in Table10.Table 10
The failure probability of constraints.
Constraints F n ≤ 0.32 25000 ≤ D W ≤ 500000 T - 0.45 D W 0.31 ≤ 0 pf MC (100) 0% 47% 11.2% PCE (20) 0% 45.2% 10.6%As shown in Table10, taking the influence of uncertainties into account, the deterministic optimal solution will violate the constraint, resulting in the failure of the solution; thus it is necessary to consider the uncertain factor in the ship design. Meanwhile, compared with the result of MC method, the error of the failure probability obtained by PCE is small; therefore, the proposed method can replace MC method to solve the failure probability.Parameter Setting of UDO. All 3 types of UDO use NSGA-II algorithm, with the population size being 40 and the number of generations being 40. Design variables, constraints, and objective functions are presented in Table 11. Models before and after dimension reduction, namely, 5-dimensional 4th-order PCE and 3-dimensinal 3rd-order PCE are both adopted for UDO and validated by MC method. Results of these 2 methods are presented in Table 12 and the iteration of the objective is shown in Figures 9, 10, and 11. Due to the randomness of NSGA-II, there are some differences in the results of the repeated search of the optimization algorithm for the same optimization problem. In this paper, the average of 5 independent times of optimization is used as the optimal value.Table 11
Design variables, uncertain variables, constraints, and objectives of UDO.
Type of UDO RDO RBDO RBRDO Design variables Design variables in Table7  Uncertain variables X : μ X = 1.7 , σ X = 0.1 μ X Y : μ Y = 0.7 , σ Y = 0.1 μ Y Z : μ Z = 0 , σ Z = 0.1 V k  Constraints Constraints in Table7 p f ( F n ≤ 0.32 ) ≤ 1 % p f ( 25000 ≤ D W ≤ 500000 ) ≤ 1 % p f ( T - 0.45 D W 0.31 ≤ 0 ) ≤ 1 %  Objectives m i n  μ ( T C ) / σ ( T C ) m i n  T C m i n  μ ( T C ) / σ ( T C )Table 12
Optimization results.
(a) Initial DO RDO PCE(5D) PCE(3D) MC(3D) T C / ( P o u n d / y r ) 9.926 8.446 8.729 8.488 8.495 M e a n _ T C / ( P o u n d / y r ) - 8.624 9.906 9.788 9.800 S T D _ T C - 4.290 4.109 3.552 3.567 L / m 195.000 188.610 181.121 178.329 178.612 B / m 32.310 31.340 32.305 29.699 29.621 D / m 20.000 15.960 15.772 15.734 15.733 T / m 10.500 11.710 11.705 11.710 11.707 C B 0.700 0.640 0.748 0.750 0.749 V k / k n 16.000 14.210 14.425 14.002 14.010 p f ( F n ≤ 0.32 )   0 0 0 p f ( 25000 ≤ D W ≤ 500000 )   41.2% 34.3% 34.3% p f ( T - 0.45 D W 0.31 ≤ 0 )   17.9% 12.1% 10.8%  (b) RBDO RBRDO PCE(5D) PCE(3D) MC(3D) PCE(5D) PCE(3D) MC(3D) T C / ( P o u n d / y r ) 8.811 8.576 8.583 8.843 8.603 8.639 M e a n _ T C / ( P o u n d / y r ) 10.128 10.077 10.091 10.227 10.042 10.019 S T D _ T C 4.237 4.128 4.147 4.229 4.088 4.077 L / m 205.357 197.358 197.873 200.666 195.025 191.701 B / m 32.310 32.309 32.253 31.420 32.299 31.786 D / m 15.719 14.521 14.529 15.784 14.518 14.316 T / m 10.769 10.859 10.854 11.460 10.831 10.704 C B 0.747 0.749 0.750 0.745 0.749 0.750 V k / k n 14.205 14.004 14.007 14.416 14.083 14.021 p f ( F n ≤ 0.32 ) 0 0 0 0 0 0 p f ( 25000 ≤ D W ≤ 500000 ) 0 0 0 0.5% 0.4% 0.1% p f ( T - 0.45 D W 0.31 ≤ 0 ) 0 0 0 0.6% 0.3% 0.1%Figure 9
Comparison of objective values in RDO process.Figure 10
Comparison of objective values in RBDO process.Figure 11
Comparison of objective values in RBRDO process.


#### 5.3.3. Analyses of Optimization Results

Some conclusions can be drawn.(1) As can be seen from Figures9, 10, and 11, for RDO, 5D PCE model starts to converge in the 900th iteration while 3D PCE model begins to converge in the 400th iteration, which indicates that the model after dimension reduction converges faster than the original model. For RBDO and RBRDO, due to the complexity of solving the failure probability, the convergence rate of the model after dimension reduction is nearly the same as the original model. However, from the total time of completing the whole optimization process, due to the reason that the collocation points of the model after dimension reduction have been reduced from 126 to 21, therefore, the time of completing 1600 iterations is much more less than that of original models, reflecting remarkable effects of dimension reduction. Besides, the difference of objective values before and after dimension reduction is little; concretely, the error of 3 cases is 1.2%, 2.7%, and 1.8%, respectively. The main reason is that TSI of A and B is relatively small, which can be regarded as uninfluential variables; therefore, the effects of fixing A and B on the optimal solution are extremely limited and the quality of the optimal solution is not greatly affected. In conclusion, the model after dimension reduction can replace the original model for UDO.(2) Table12 also gives the optimal solutions after dimension reduction (3D) of RDO, RBDO, and RBRDO, solved by PCE and MC method. Optimal cases obtained by these 2 methods are almost the same, which validates the feasibility of PCE again. Moreover, the advantage of PCE method can be seen from Figures 9, 10, and 11; due to the fact that PCE can maintain the accuracy with fewer sample points, for RDO, the objective of PCE has decreased greatly in the early stage of the iteration and tends to converge at the 400th iteration. The objective of PCE completes the entire 1600 iterations in 400 min while the objective of MC method completes the entire 1600 iterations in 1800 min. Similarly, for RBDO, the objective of PCE completes the entire 1600 iterations in 1800 min while the objective of MC method completes the entire 1600 iterations in 2600 min. For RBRDO, the difference is more obvious that the objective of PCE completes the entire 1600 iterations in 3600 min while the objective of MC method completes the entire 1600 iterations in 5600 min. The time of PCE is significantly lower than that of MC method for all 3 types of UDO, thereby showing obvious advantages in engineering applications.(3) Taking the result of PCE as an example, compared with the initial design,TC of the 3 uncertain optimal cases is reduced by 14.5% (RDO), 13.6% (RBDO), and 13.2% (RBRDO), which are slightly lower than that of DO, 14.9%. However, for the 3 uncertain optimal cases, the influence of the uncertainty is taken into account during the optimization process; thus the standard deviation is reduced greatly compared with that of DO. More specifically, for RDO, its optimal case has lower mean and standard deviation than those of RBDO, which means the perturbation of the objective caused by the uncertainty is much smaller than that of the latter, leading to a more robust design. For RBDO, its failure probability of constraints is all 0, which is lower than that of RDO, thereby leading to a design with a higher reliability. For RBRDO, it combines the advantages of RDO and RBDO, the mean and the standard deviation are reduced, and the failure probability of constraints is less than 0.5%, which shows the fact that the optimal case of RBRDO can guarantee the reliability demand and reduce the perturbation caused by the uncertainty at the same time, thus ensuring both the reliability and the robustness of the case.


## 5.1. Description of Optimization Design

The proposed method is applied to the ship UDO. The optimization model of the bulk carrier proposed by Hannapel S and Vlahopoulos N is considered. Details of this model can be found in [32]. Deterministic optimization (DO) and uncertainty optimization are both applied to this model.


## 5.2. Specific Process of Bulk Carrier Uncertainty-Based Design Optimization

The specific process is shown in Figure6. The steps of this process are as follows.Figure 6
Flow chart of ship uncertainty-based design optimization.Step 1.
Preparation work: determine the optimal PCE order of original model, then perform the sensitivity analysis on uncertain variables and rank them according to their degree of influence on the objective and the constraint, fix the uncertain variable with less influence at the mean value, and reduce the dimension of the model. Redetermine the optimal PCE order for the simplified model.Step 2.
Design variables and their ranges are determined.Step 3.
A set of design variables, in other words, a bulk carrier scheme, is selected and transferred to the inner layer, the uncertainty analysis module.Step 4.
According to the expansion order, the distribution type, and the number of uncertain variables, which are determined in advance, IPCM is used to select collocation points for uncertain variables.Step 5.
PCE is used to perform the expansion of the objective and the constraint.Step 6.
According to the design requirement, select one of 3 types of UDO. Then the mean and the standard deviation of the objective and the failure probability of the constraint are obtained and transferred to the outer layer, optimization module.Step 7.
The optimizer judges whether the entire optimization process meets requirements of the iteration number. If the process reaches the iteration number, the optimization process ends. If not, step 3~step 6 are implemented again.


## 5.3. Optimization Analyses

### 5.3.1. Deterministic Optimization

In order to give a baseline for the comparison for later results of uncertainty optimization, DO is done first. DO uses a multi-island genetic algorithm, with the subpopulation size being 10, the number of islands being 10, and the number of generations being 50. Design variables, constraints, and objective functions are presented in Table7 and the comparison between results of DO and initial design [32] is listed in Table 12.Table 7
Bulk carrier model.
Design variables L e n g t h , L / m 150 ≤ L ≤ 274.32 B r e a d t h , B / m 20 ≤ B ≤ 32.31 D e p t h , D / m 13 ≤ D ≤ 25 D r a f t , T / m 10 ≤ D ≤ 11.71 B l o c k   c o e f f i c i e n t , C B 0.63 ≤ C B ≤ 0.75 D e s i g n   s p e e d , V k 14 ≤ V k ≤ 20  Objective function  T r a n s p o r t a t i o n   c o s t , T C / ( P o u n d / y r ) T C = C A / A C  constraints  L / B ≥ 6 L / D ≤ 15 L / T ≤ 19 F n ≤ 0.32 25000 ≤ D W ≤ 500000 T - 0.45 D W 0.31 ≤ 0 T - 0.7 D - 0.7 ≤ 0 0.07 B - 0.35 T - ( 0.085 C B - 0.002 ) B 2 T C B + 1 + 0.25 D ≤ 0


### 5.3.2. Uncertainty-Based Optimization

Selection of Uncertain Variables. In the real travelling process, due to the unavoidable disturbance such as the wind and the wave, the ship resistance in waves will increase with more violent ship motions; therefore the service/operational speed of ship will decrease compared with the hydrostatic condition. Therefore, the speed of a sailing ship is not a fixed value, but an uncertain variable that perturbs around the design speed. By monitoring the ship speed in a period of time [6, 13], we obviously find that the speed perturbation is a universal existence. Although the perturbation uncertainty magnitudes are relatively small in most cases, their continuous calculation and coupling effects with other parameters will cause the system response to a large deviation in the optimization design process. Therefore, the ship speed is selected as an uncertain variable in (32)Vactual=Vk+ZwhereVactual is the actual speed, Vk is the design speed, and Z is the speed perturbation, which follows the normal distribution with the mean of 0 and the standard deviation of 0.1Vk [6].Then the uncertainty in parameters is introduced. The steel weight of the bulk carrier is given by(33)WS=0.034L1.7B0.7D0.4CB0.5Equation (33) is obtained by the statistical regression using the real ship data. Considering the limited number of samples and the uncontrollable error of using the unknown regression method, the regression exponent cannot accurately reflect the actual engineering situation; therefore, the uncertainty is applied to the exponent [32], which are assumed to be random parameters, namely, X (the normal distribution with the mean of 1.7 and the standard deviation of 0.17) and Y (the normal distribution with the mean of 0.7 and the standard deviation of 0.07), A (the normal distribution with the mean of 0.4 and the standard deviation of 0.04), and B (the normal distribution with the mean of 0.5 and the standard deviation of 0.05), as presented in (34)WS=0.034LXBYDACBBDetermination of the Optimal PCE Order. Considering the influence of the uncertainty on the objective and the constraint for the DO optimal case, the objective and the constraint are computed using PCE, whose expansion order is d=2,3,4 with the number of their corresponding expansion terms P are 21, 56, and 126, respectively. IPCM is used to select collocation points, for the optimization objective, and the constraint, the value of undetermined coefficients of order d=2,3,4 is shown in Figure 7, where the abscissa is the serial number of the coefficient and the ordinate is the coefficient value.Figure 7
Coefficient values of the objective and the constraint with different expansion orders.It can be seen from Figure10 that, for T-0.45DW0.31, undetermined coefficient values after the number 72 are almost all 0; thus the 4th order expansion can basically meet the requirement of the accuracy. For Fn, the coefficient values of d=2,3,4 are nearly the same, thus the 2nd-order expansion is optimal for Fn. For DW, the coefficient values after the number 56 (the 4th terms) are almost all 0, which means the 4th terms have no contribution to DW; therefore, 3rd-order expansion is selected. For TC, coefficient values have been close to 0 since the number 107; therefore, the 4th-order can basically meet requirement of the accuracy. In summary, the optimal PCE order is the 4th-order.Sensitivity Analysis. Then, for the objective and the constraint, sensitivity analysis is performed by using PCE method and traditional Sobol’ method. The sensitivity indices, namely, MSI and TSI of the objective and the constraint to X, Y, Z, A, and B, are shown in Table 8 and Figure 8.Table 8
Sensitivity indices of the objective and constraints.
T - 0.45 D W 0.31 D W F n T C PCE Sobol’ PCE Sobol’ PCE Sobol PCE Sobol’ M S I X 0.8692 0.8857 0.8842 0.8302 0 0 0.8507 0.8925 M S I Y 0.0615 0.0637 0.0595 0.0582 0 0 0.0649 0.0712 M S I Z 0.0002 0.0012 0.0002 0.0004 0.9999 0.9249 0.0023 0.0021 M S I A 0.0136 0.0153 0.0128 0.0139 0 0 0.0148 0.151 M S I B 0.0006 0.0005 0.0005 0.0005 0 0 0.0007 0.0006 T S I X 0.9238 0.9159 0.9267 0.9451 0 0 0.9169 0.9237 T S I Y 0.1060 0.1155 0.0945 0.1151 0 0 0.1180 0.1244 T S I Z 0.0002 0.0003 0.0002 0.0002 0.9999 0.9249 0.0027 0.0031 T S I A 0.0241 0.0230 0.0208 0.0314 0 0 0.0279 0.0293 T S I B 0.0011 0.0020 0.0009 0.0012 0 0 0.0013 0.0021Figure 8
Sensitivity indices of the objective and constraints.According toMSI and TSI, 5 uncertain variables are ranked in order. For T-0.45DW0.31, DW and TC, X, and Y have dominated the role and A has only a little effect, while the influence of other variables can be neglected. On the contrary, for Fn, Z appears to be the most important in the decomposition of the variance of the response while the effect of other 4 variables has no need to consider, which is consistent with the formula of Fn. In conclusion, whether from MSI or TSI, X, Y, and Z have main effects on the objective and the constraint. Therefore, A and B are fixed at their mean values, while X, Y, and Z still satisfy their original distributions. The simplified model is used to replace the original model to perform optimization.For this simplified model, the process of determining the optimal PCE order is repeated, and 3rd-order expansion is optimal for the model. The mean and the standard deviation of 3rd-order expansion are compared with results of MC method, as can be seen in Table9.Table 9
The mean and the standard deviation of the objective and the constraint.
the objective /the constraint T C F n ≤ 0.32 25000 ≤ D W ≤ 500000 T - 0.45 D W 0.31 ≤ 0 Mean MC (100) 10.06 0.1699 3.408∗104 0.309 PCE (20) 9.92 0.1699 3.479∗104 0.304 Error 1.39 % 0 % 1.95 % 1.64 %  Standard deviation MC (100) 4.21 0.017 8009 0.988 PCE (20) 4.29 0.017 7968 1.021 Error 1.90 % 0 % 0.52 % 3.18 %From results in Table9, taking the result of MC method as the exact value, the error of the mean and the standard deviation solved by PCE is within 2% and 3.5%, respectively. Meanwhile required sample points of PCE (20 sample points) are far less than those of MC method (100 sample points). Therefore, the advantage of PCE in computational efforts is obvious; specifically, required sample points reduce greatly while the calculation accuracy is not significantly reduced.Then, according to the solved first 4 order moments of the constraints, maximum entropy method (MEM) is adopted to obtain PDF of constraints and then the failure probability is solved. Compared with MC method, the results are listed in Table10.Table 10
The failure probability of constraints.
Constraints F n ≤ 0.32 25000 ≤ D W ≤ 500000 T - 0.45 D W 0.31 ≤ 0 pf MC (100) 0% 47% 11.2% PCE (20) 0% 45.2% 10.6%As shown in Table10, taking the influence of uncertainties into account, the deterministic optimal solution will violate the constraint, resulting in the failure of the solution; thus it is necessary to consider the uncertain factor in the ship design. Meanwhile, compared with the result of MC method, the error of the failure probability obtained by PCE is small; therefore, the proposed method can replace MC method to solve the failure probability.Parameter Setting of UDO. All 3 types of UDO use NSGA-II algorithm, with the population size being 40 and the number of generations being 40. Design variables, constraints, and objective functions are presented in Table 11. Models before and after dimension reduction, namely, 5-dimensional 4th-order PCE and 3-dimensinal 3rd-order PCE are both adopted for UDO and validated by MC method. Results of these 2 methods are presented in Table 12 and the iteration of the objective is shown in Figures 9, 10, and 11. Due to the randomness of NSGA-II, there are some differences in the results of the repeated search of the optimization algorithm for the same optimization problem. In this paper, the average of 5 independent times of optimization is used as the optimal value.Table 11
Design variables, uncertain variables, constraints, and objectives of UDO.
Type of UDO RDO RBDO RBRDO Design variables Design variables in Table7  Uncertain variables X : μ X = 1.7 , σ X = 0.1 μ X Y : μ Y = 0.7 , σ Y = 0.1 μ Y Z : μ Z = 0 , σ Z = 0.1 V k  Constraints Constraints in Table7 p f ( F n ≤ 0.32 ) ≤ 1 % p f ( 25000 ≤ D W ≤ 500000 ) ≤ 1 % p f ( T - 0.45 D W 0.31 ≤ 0 ) ≤ 1 %  Objectives m i n  μ ( T C ) / σ ( T C ) m i n  T C m i n  μ ( T C ) / σ ( T C )Table 12
Optimization results.
(a) Initial DO RDO PCE(5D) PCE(3D) MC(3D) T C / ( P o u n d / y r ) 9.926 8.446 8.729 8.488 8.495 M e a n _ T C / ( P o u n d / y r ) - 8.624 9.906 9.788 9.800 S T D _ T C - 4.290 4.109 3.552 3.567 L / m 195.000 188.610 181.121 178.329 178.612 B / m 32.310 31.340 32.305 29.699 29.621 D / m 20.000 15.960 15.772 15.734 15.733 T / m 10.500 11.710 11.705 11.710 11.707 C B 0.700 0.640 0.748 0.750 0.749 V k / k n 16.000 14.210 14.425 14.002 14.010 p f ( F n ≤ 0.32 )   0 0 0 p f ( 25000 ≤ D W ≤ 500000 )   41.2% 34.3% 34.3% p f ( T - 0.45 D W 0.31 ≤ 0 )   17.9% 12.1% 10.8%  (b) RBDO RBRDO PCE(5D) PCE(3D) MC(3D) PCE(5D) PCE(3D) MC(3D) T C / ( P o u n d / y r ) 8.811 8.576 8.583 8.843 8.603 8.639 M e a n _ T C / ( P o u n d / y r ) 10.128 10.077 10.091 10.227 10.042 10.019 S T D _ T C 4.237 4.128 4.147 4.229 4.088 4.077 L / m 205.357 197.358 197.873 200.666 195.025 191.701 B / m 32.310 32.309 32.253 31.420 32.299 31.786 D / m 15.719 14.521 14.529 15.784 14.518 14.316 T / m 10.769 10.859 10.854 11.460 10.831 10.704 C B 0.747 0.749 0.750 0.745 0.749 0.750 V k / k n 14.205 14.004 14.007 14.416 14.083 14.021 p f ( F n ≤ 0.32 ) 0 0 0 0 0 0 p f ( 25000 ≤ D W ≤ 500000 ) 0 0 0 0.5% 0.4% 0.1% p f ( T - 0.45 D W 0.31 ≤ 0 ) 0 0 0 0.6% 0.3% 0.1%Figure 9
Comparison of objective values in RDO process.Figure 10
Comparison of objective values in RBDO process.Figure 11
Comparison of objective values in RBRDO process.


### 5.3.3. Analyses of Optimization Results

Some conclusions can be drawn.(1) As can be seen from Figures9, 10, and 11, for RDO, 5D PCE model starts to converge in the 900th iteration while 3D PCE model begins to converge in the 400th iteration, which indicates that the model after dimension reduction converges faster than the original model. For RBDO and RBRDO, due to the complexity of solving the failure probability, the convergence rate of the model after dimension reduction is nearly the same as the original model. However, from the total time of completing the whole optimization process, due to the reason that the collocation points of the model after dimension reduction have been reduced from 126 to 21, therefore, the time of completing 1600 iterations is much more less than that of original models, reflecting remarkable effects of dimension reduction. Besides, the difference of objective values before and after dimension reduction is little; concretely, the error of 3 cases is 1.2%, 2.7%, and 1.8%, respectively. The main reason is that TSI of A and B is relatively small, which can be regarded as uninfluential variables; therefore, the effects of fixing A and B on the optimal solution are extremely limited and the quality of the optimal solution is not greatly affected. In conclusion, the model after dimension reduction can replace the original model for UDO.(2) Table12 also gives the optimal solutions after dimension reduction (3D) of RDO, RBDO, and RBRDO, solved by PCE and MC method. Optimal cases obtained by these 2 methods are almost the same, which validates the feasibility of PCE again. Moreover, the advantage of PCE method can be seen from Figures 9, 10, and 11; due to the fact that PCE can maintain the accuracy with fewer sample points, for RDO, the objective of PCE has decreased greatly in the early stage of the iteration and tends to converge at the 400th iteration. The objective of PCE completes the entire 1600 iterations in 400 min while the objective of MC method completes the entire 1600 iterations in 1800 min. Similarly, for RBDO, the objective of PCE completes the entire 1600 iterations in 1800 min while the objective of MC method completes the entire 1600 iterations in 2600 min. For RBRDO, the difference is more obvious that the objective of PCE completes the entire 1600 iterations in 3600 min while the objective of MC method completes the entire 1600 iterations in 5600 min. The time of PCE is significantly lower than that of MC method for all 3 types of UDO, thereby showing obvious advantages in engineering applications.(3) Taking the result of PCE as an example, compared with the initial design,TC of the 3 uncertain optimal cases is reduced by 14.5% (RDO), 13.6% (RBDO), and 13.2% (RBRDO), which are slightly lower than that of DO, 14.9%. However, for the 3 uncertain optimal cases, the influence of the uncertainty is taken into account during the optimization process; thus the standard deviation is reduced greatly compared with that of DO. More specifically, for RDO, its optimal case has lower mean and standard deviation than those of RBDO, which means the perturbation of the objective caused by the uncertainty is much smaller than that of the latter, leading to a more robust design. For RBDO, its failure probability of constraints is all 0, which is lower than that of RDO, thereby leading to a design with a higher reliability. For RBRDO, it combines the advantages of RDO and RBDO, the mean and the standard deviation are reduced, and the failure probability of constraints is less than 0.5%, which shows the fact that the optimal case of RBRDO can guarantee the reliability demand and reduce the perturbation caused by the uncertainty at the same time, thus ensuring both the reliability and the robustness of the case.


## 5.3.1. Deterministic Optimization

In order to give a baseline for the comparison for later results of uncertainty optimization, DO is done first. DO uses a multi-island genetic algorithm, with the subpopulation size being 10, the number of islands being 10, and the number of generations being 50. Design variables, constraints, and objective functions are presented in Table7 and the comparison between results of DO and initial design [32] is listed in Table 12.Table 7
Bulk carrier model.
Design variables L e n g t h , L / m 150 ≤ L ≤ 274.32 B r e a d t h , B / m 20 ≤ B ≤ 32.31 D e p t h , D / m 13 ≤ D ≤ 25 D r a f t , T / m 10 ≤ D ≤ 11.71 B l o c k   c o e f f i c i e n t , C B 0.63 ≤ C B ≤ 0.75 D e s i g n   s p e e d , V k 14 ≤ V k ≤ 20  Objective function  T r a n s p o r t a t i o n   c o s t , T C / ( P o u n d / y r ) T C = C A / A C  constraints  L / B ≥ 6 L / D ≤ 15 L / T ≤ 19 F n ≤ 0.32 25000 ≤ D W ≤ 500000 T - 0.45 D W 0.31 ≤ 0 T - 0.7 D - 0.7 ≤ 0 0.07 B - 0.35 T - ( 0.085 C B - 0.002 ) B 2 T C B + 1 + 0.25 D ≤ 0


## 5.3.2. Uncertainty-Based Optimization

Selection of Uncertain Variables. In the real travelling process, due to the unavoidable disturbance such as the wind and the wave, the ship resistance in waves will increase with more violent ship motions; therefore the service/operational speed of ship will decrease compared with the hydrostatic condition. Therefore, the speed of a sailing ship is not a fixed value, but an uncertain variable that perturbs around the design speed. By monitoring the ship speed in a period of time [6, 13], we obviously find that the speed perturbation is a universal existence. Although the perturbation uncertainty magnitudes are relatively small in most cases, their continuous calculation and coupling effects with other parameters will cause the system response to a large deviation in the optimization design process. Therefore, the ship speed is selected as an uncertain variable in (32)Vactual=Vk+ZwhereVactual is the actual speed, Vk is the design speed, and Z is the speed perturbation, which follows the normal distribution with the mean of 0 and the standard deviation of 0.1Vk [6].Then the uncertainty in parameters is introduced. The steel weight of the bulk carrier is given by(33)WS=0.034L1.7B0.7D0.4CB0.5Equation (33) is obtained by the statistical regression using the real ship data. Considering the limited number of samples and the uncontrollable error of using the unknown regression method, the regression exponent cannot accurately reflect the actual engineering situation; therefore, the uncertainty is applied to the exponent [32], which are assumed to be random parameters, namely, X (the normal distribution with the mean of 1.7 and the standard deviation of 0.17) and Y (the normal distribution with the mean of 0.7 and the standard deviation of 0.07), A (the normal distribution with the mean of 0.4 and the standard deviation of 0.04), and B (the normal distribution with the mean of 0.5 and the standard deviation of 0.05), as presented in (34)WS=0.034LXBYDACBBDetermination of the Optimal PCE Order. Considering the influence of the uncertainty on the objective and the constraint for the DO optimal case, the objective and the constraint are computed using PCE, whose expansion order is d=2,3,4 with the number of their corresponding expansion terms P are 21, 56, and 126, respectively. IPCM is used to select collocation points, for the optimization objective, and the constraint, the value of undetermined coefficients of order d=2,3,4 is shown in Figure 7, where the abscissa is the serial number of the coefficient and the ordinate is the coefficient value.Figure 7
Coefficient values of the objective and the constraint with different expansion orders.It can be seen from Figure10 that, for T-0.45DW0.31, undetermined coefficient values after the number 72 are almost all 0; thus the 4th order expansion can basically meet the requirement of the accuracy. For Fn, the coefficient values of d=2,3,4 are nearly the same, thus the 2nd-order expansion is optimal for Fn. For DW, the coefficient values after the number 56 (the 4th terms) are almost all 0, which means the 4th terms have no contribution to DW; therefore, 3rd-order expansion is selected. For TC, coefficient values have been close to 0 since the number 107; therefore, the 4th-order can basically meet requirement of the accuracy. In summary, the optimal PCE order is the 4th-order.Sensitivity Analysis. Then, for the objective and the constraint, sensitivity analysis is performed by using PCE method and traditional Sobol’ method. The sensitivity indices, namely, MSI and TSI of the objective and the constraint to X, Y, Z, A, and B, are shown in Table 8 and Figure 8.Table 8
Sensitivity indices of the objective and constraints.
T - 0.45 D W 0.31 D W F n T C PCE Sobol’ PCE Sobol’ PCE Sobol PCE Sobol’ M S I X 0.8692 0.8857 0.8842 0.8302 0 0 0.8507 0.8925 M S I Y 0.0615 0.0637 0.0595 0.0582 0 0 0.0649 0.0712 M S I Z 0.0002 0.0012 0.0002 0.0004 0.9999 0.9249 0.0023 0.0021 M S I A 0.0136 0.0153 0.0128 0.0139 0 0 0.0148 0.151 M S I B 0.0006 0.0005 0.0005 0.0005 0 0 0.0007 0.0006 T S I X 0.9238 0.9159 0.9267 0.9451 0 0 0.9169 0.9237 T S I Y 0.1060 0.1155 0.0945 0.1151 0 0 0.1180 0.1244 T S I Z 0.0002 0.0003 0.0002 0.0002 0.9999 0.9249 0.0027 0.0031 T S I A 0.0241 0.0230 0.0208 0.0314 0 0 0.0279 0.0293 T S I B 0.0011 0.0020 0.0009 0.0012 0 0 0.0013 0.0021Figure 8
Sensitivity indices of the objective and constraints.According toMSI and TSI, 5 uncertain variables are ranked in order. For T-0.45DW0.31, DW and TC, X, and Y have dominated the role and A has only a little effect, while the influence of other variables can be neglected. On the contrary, for Fn, Z appears to be the most important in the decomposition of the variance of the response while the effect of other 4 variables has no need to consider, which is consistent with the formula of Fn. In conclusion, whether from MSI or TSI, X, Y, and Z have main effects on the objective and the constraint. Therefore, A and B are fixed at their mean values, while X, Y, and Z still satisfy their original distributions. The simplified model is used to replace the original model to perform optimization.For this simplified model, the process of determining the optimal PCE order is repeated, and 3rd-order expansion is optimal for the model. The mean and the standard deviation of 3rd-order expansion are compared with results of MC method, as can be seen in Table9.Table 9
The mean and the standard deviation of the objective and the constraint.
the objective /the constraint T C F n ≤ 0.32 25000 ≤ D W ≤ 500000 T - 0.45 D W 0.31 ≤ 0 Mean MC (100) 10.06 0.1699 3.408∗104 0.309 PCE (20) 9.92 0.1699 3.479∗104 0.304 Error 1.39 % 0 % 1.95 % 1.64 %  Standard deviation MC (100) 4.21 0.017 8009 0.988 PCE (20) 4.29 0.017 7968 1.021 Error 1.90 % 0 % 0.52 % 3.18 %From results in Table9, taking the result of MC method as the exact value, the error of the mean and the standard deviation solved by PCE is within 2% and 3.5%, respectively. Meanwhile required sample points of PCE (20 sample points) are far less than those of MC method (100 sample points). Therefore, the advantage of PCE in computational efforts is obvious; specifically, required sample points reduce greatly while the calculation accuracy is not significantly reduced.Then, according to the solved first 4 order moments of the constraints, maximum entropy method (MEM) is adopted to obtain PDF of constraints and then the failure probability is solved. Compared with MC method, the results are listed in Table10.Table 10
The failure probability of constraints.
Constraints F n ≤ 0.32 25000 ≤ D W ≤ 500000 T - 0.45 D W 0.31 ≤ 0 pf MC (100) 0% 47% 11.2% PCE (20) 0% 45.2% 10.6%As shown in Table10, taking the influence of uncertainties into account, the deterministic optimal solution will violate the constraint, resulting in the failure of the solution; thus it is necessary to consider the uncertain factor in the ship design. Meanwhile, compared with the result of MC method, the error of the failure probability obtained by PCE is small; therefore, the proposed method can replace MC method to solve the failure probability.Parameter Setting of UDO. All 3 types of UDO use NSGA-II algorithm, with the population size being 40 and the number of generations being 40. Design variables, constraints, and objective functions are presented in Table 11. Models before and after dimension reduction, namely, 5-dimensional 4th-order PCE and 3-dimensinal 3rd-order PCE are both adopted for UDO and validated by MC method. Results of these 2 methods are presented in Table 12 and the iteration of the objective is shown in Figures 9, 10, and 11. Due to the randomness of NSGA-II, there are some differences in the results of the repeated search of the optimization algorithm for the same optimization problem. In this paper, the average of 5 independent times of optimization is used as the optimal value.Table 11
Design variables, uncertain variables, constraints, and objectives of UDO.
Type of UDO RDO RBDO RBRDO Design variables Design variables in Table7  Uncertain variables X : μ X = 1.7 , σ X = 0.1 μ X Y : μ Y = 0.7 , σ Y = 0.1 μ Y Z : μ Z = 0 , σ Z = 0.1 V k  Constraints Constraints in Table7 p f ( F n ≤ 0.32 ) ≤ 1 % p f ( 25000 ≤ D W ≤ 500000 ) ≤ 1 % p f ( T - 0.45 D W 0.31 ≤ 0 ) ≤ 1 %  Objectives m i n  μ ( T C ) / σ ( T C ) m i n  T C m i n  μ ( T C ) / σ ( T C )Table 12
Optimization results.
(a) Initial DO RDO PCE(5D) PCE(3D) MC(3D) T C / ( P o u n d / y r ) 9.926 8.446 8.729 8.488 8.495 M e a n _ T C / ( P o u n d / y r ) - 8.624 9.906 9.788 9.800 S T D _ T C - 4.290 4.109 3.552 3.567 L / m 195.000 188.610 181.121 178.329 178.612 B / m 32.310 31.340 32.305 29.699 29.621 D / m 20.000 15.960 15.772 15.734 15.733 T / m 10.500 11.710 11.705 11.710 11.707 C B 0.700 0.640 0.748 0.750 0.749 V k / k n 16.000 14.210 14.425 14.002 14.010 p f ( F n ≤ 0.32 )   0 0 0 p f ( 25000 ≤ D W ≤ 500000 )   41.2% 34.3% 34.3% p f ( T - 0.45 D W 0.31 ≤ 0 )   17.9% 12.1% 10.8%  (b) RBDO RBRDO PCE(5D) PCE(3D) MC(3D) PCE(5D) PCE(3D) MC(3D) T C / ( P o u n d / y r ) 8.811 8.576 8.583 8.843 8.603 8.639 M e a n _ T C / ( P o u n d / y r ) 10.128 10.077 10.091 10.227 10.042 10.019 S T D _ T C 4.237 4.128 4.147 4.229 4.088 4.077 L / m 205.357 197.358 197.873 200.666 195.025 191.701 B / m 32.310 32.309 32.253 31.420 32.299 31.786 D / m 15.719 14.521 14.529 15.784 14.518 14.316 T / m 10.769 10.859 10.854 11.460 10.831 10.704 C B 0.747 0.749 0.750 0.745 0.749 0.750 V k / k n 14.205 14.004 14.007 14.416 14.083 14.021 p f ( F n ≤ 0.32 ) 0 0 0 0 0 0 p f ( 25000 ≤ D W ≤ 500000 ) 0 0 0 0.5% 0.4% 0.1% p f ( T - 0.45 D W 0.31 ≤ 0 ) 0 0 0 0.6% 0.3% 0.1%Figure 9
Comparison of objective values in RDO process.Figure 10
Comparison of objective values in RBDO process.Figure 11
Comparison of objective values in RBRDO process.


## 5.3.3. Analyses of Optimization Results

Some conclusions can be drawn.(1) As can be seen from Figures9, 10, and 11, for RDO, 5D PCE model starts to converge in the 900th iteration while 3D PCE model begins to converge in the 400th iteration, which indicates that the model after dimension reduction converges faster than the original model. For RBDO and RBRDO, due to the complexity of solving the failure probability, the convergence rate of the model after dimension reduction is nearly the same as the original model. However, from the total time of completing the whole optimization process, due to the reason that the collocation points of the model after dimension reduction have been reduced from 126 to 21, therefore, the time of completing 1600 iterations is much more less than that of original models, reflecting remarkable effects of dimension reduction. Besides, the difference of objective values before and after dimension reduction is little; concretely, the error of 3 cases is 1.2%, 2.7%, and 1.8%, respectively. The main reason is that TSI of A and B is relatively small, which can be regarded as uninfluential variables; therefore, the effects of fixing A and B on the optimal solution are extremely limited and the quality of the optimal solution is not greatly affected. In conclusion, the model after dimension reduction can replace the original model for UDO.(2) Table12 also gives the optimal solutions after dimension reduction (3D) of RDO, RBDO, and RBRDO, solved by PCE and MC method. Optimal cases obtained by these 2 methods are almost the same, which validates the feasibility of PCE again. Moreover, the advantage of PCE method can be seen from Figures 9, 10, and 11; due to the fact that PCE can maintain the accuracy with fewer sample points, for RDO, the objective of PCE has decreased greatly in the early stage of the iteration and tends to converge at the 400th iteration. The objective of PCE completes the entire 1600 iterations in 400 min while the objective of MC method completes the entire 1600 iterations in 1800 min. Similarly, for RBDO, the objective of PCE completes the entire 1600 iterations in 1800 min while the objective of MC method completes the entire 1600 iterations in 2600 min. For RBRDO, the difference is more obvious that the objective of PCE completes the entire 1600 iterations in 3600 min while the objective of MC method completes the entire 1600 iterations in 5600 min. The time of PCE is significantly lower than that of MC method for all 3 types of UDO, thereby showing obvious advantages in engineering applications.(3) Taking the result of PCE as an example, compared with the initial design,TC of the 3 uncertain optimal cases is reduced by 14.5% (RDO), 13.6% (RBDO), and 13.2% (RBRDO), which are slightly lower than that of DO, 14.9%. However, for the 3 uncertain optimal cases, the influence of the uncertainty is taken into account during the optimization process; thus the standard deviation is reduced greatly compared with that of DO. More specifically, for RDO, its optimal case has lower mean and standard deviation than those of RBDO, which means the perturbation of the objective caused by the uncertainty is much smaller than that of the latter, leading to a more robust design. For RBDO, its failure probability of constraints is all 0, which is lower than that of RDO, thereby leading to a design with a higher reliability. For RBRDO, it combines the advantages of RDO and RBDO, the mean and the standard deviation are reduced, and the failure probability of constraints is less than 0.5%, which shows the fact that the optimal case of RBRDO can guarantee the reliability demand and reduce the perturbation caused by the uncertainty at the same time, thus ensuring both the reliability and the robustness of the case.


## 6. Conclusions

This paper proposed a sensitivity analysis based on polynomial chaos expansions and its application in ship UDO. Combined with the multiobjective genetic algorithm, an efficient UDO system is constructed and applied to the ship design, and some conclusions can be drawn:(1) PCE-based Sobol’ indices are used to perform SA for the model and uncertain variables are ranked based on the indices to achieve the dimensionality reduction. Compared with MC method, the proposed method has proven the efficiency and accuracy in SA. Specifically speaking, the computation of Sobol’ indices after a proper expansion is precise enough while the computation cost is transferred to the obtention of the PCE coefficients, the subsequent postprocessing being almost costless.(2) When using PCE, IPCM can give the optimal number of collocation points by comparing the rank of the coefficient matrix with the number of collocation points, which reduces collocation points greatly and overcomes the shortcoming, large computational complexity of traditional statistical methods.(3) In the bulk carrier uncertainty-based design optimization, the perturbation of the objective and the constraint under the influence of multiple random variables are considered. By analyzing optimization results, first, it is necessary to take the influence of uncertainties into account where the deterministic optimal solution will violate the constraint, resulting in the failure of the solution. Meanwhile, in the optimization process, on the premise of maintaining the accuracy, the time of models after dimension reduction is significantly lower than that of original models; at the same time, the time of PCE is significantly lower than that of commonly used MC method, which shows more advantages in engineering applications.However, the paper models all uncertain variables by using only probabilistic method and only discussed the application of the proposed method in the ship preliminary design optimization. Future work will use the probabilistic and nonprobabilistic methods to model uncertain variables according to their different characteristics and will not be limited to the ship conceptual design and will be carried out to guide the hull form design.



---

*Source: 7498526-2019-01-23.xml*