Manufacturing cost and reliability‐based shape optimization of plate structures

A novel methodology is presented for the manufacturing cost and reliability‐based optimization of plate structures with the boundary element method (BEM), with the aim of improving the accuracy, robustness, and efficiency of the optimization of aircraft structures. The BEM plate formulations with respect to plate thickness are derived for the first time, and used as part of an implicit differentiation method (IDM), enabling the full shape optimization of plate structures with the BEM. These implicit derivatives are validated against derivatives obtained from the finite difference method (FDM) and from an analytical solution. Results indicate that the IDM is more robust than the FDM and in excellent agreement with the analytical solution, and more accurate than the FDM for most of the step‐sizes investigated. To demonstrate the full shape optimization of plates with the newly developed IDM, a numerical example involving reliability‐based design optimization and manufacturing cost optimization is presented for a plate structure. Results show that the newly developed IDM is more efficient than the FDM when performing this optimization.


INTRODUCTION
Structural mass is a very important consideration for aeronautical engineers, as small increases in the mass of an aircraft can lead to significant increases in fuel consumption and carbon emissions. Alongside mass, manufacturing costs are also a significant consideration for aeronautical engineers. Reductions in manufacturing costs can enable aircraft manufacturers to divert more resources towards improving safety and reducing the environmental impact of their aircraft. Reductions in costs could also be passed down to consumers via reduced ticket prices. Reliability-based design optimization (RBDO) can often provide designs for structures that are optimized in terms of mass and reliability, but not necessarily in terms of manufacturing costs, especially for complex structures. This work aims to develop a methodology that couples the RBDO of aircraft structures with a manufacturing cost estimation approach.
Optimizing a structure for manufacturing cost and reliability involves shape optimization-optimizing the shape of a structure to minimize manufacturing cost and to maximize reliability. The boundary element method (BEM) can be a very effective tool for shape optimization. The BEM only discretizes the outer boundary of a structure, allowing the outer geometry to be varied without requiring the entire structure to be resmeshed, saving both time and computational resources. An added benefit of this feature is that it allows the sensitivities of the responses of a BEM model, with respect to changes in the outer geometry, to be calculated in a very straightforward and computationally efficient manner, making the BEM a very effective tool for the shape optimization of engineering structures.
Previous work on the topic of shape optimization with the BEM by the research community has mostly considered 2D structures, [1][2][3][4][5][6][7][8][9][10][11] and to a lesser-degree 3D structures. [12][13][14][15][16][17][18] Very few previous works by the research community on the topic of shape optimization with the BEM have considered plate or shell structures. Where possible, modeling a structure as a plate or shell structure can have significant benefits over modeling the structure in 2D or 3D. This is because many real-life engineering structures cannot be accurately modeled as 2D structures, and depending on the engineering structure being modeled, it can often be more efficient to model a structure as a plate or a shell structure, rather than as a 3D structure, without a noticeable loss in modeling accuracy. An example of this could be a curved stiffened panel from an aircraft fuselage, or a wingbox from an aircraft's wing. To the authors' knowledge, there is only one previous work by the research community on the topic of shape optimization with the BEM with plate or shell structures; Babouskos et al. 19 optimized the thickness distribution in a thin plate to regulate the dynamic response of the plate. The thickness distribution of the plate was approximated via a surrogate model, and the optimal thickness distribution in the plate was calculated using derivatives of this surrogate model. The derivatives of the BEM formulations for plates were not calculated, and the remaining geometry of the plate was not optimized. This current work aims to develop an implicit differentiation method (IDM), which uses the implicit/direct derivatives of the BEM formulations for plates, to enable the full shape optimization of plate structures with the BEM, involving all of a plate's geometry. By using the implicit/direct derivatives of the BEM formulations, shape optimization can be conducted in a much more accurate and efficient manner than with other methods.
The first steps towards developing an IDM for plate or shell structures with the BEM were conducted by Morse et al. in Reference 20, in which the exact/implicit derivatives of the BEM plate formulations with respect to geometrical variables were derived for the first time. These exact derivatives were only applicable to geometrical variables that influence the nodal coordinates of the BEM plate mesh, such as plate length or width, and so they were not applicable to geometrical variables that do not influence nodal coordinates, such as plate thickness. Therefore, the full shape optimization of plate structures, involving all geometrical variables, was not possible. This current work aims to build upon the work presented in Reference 20 and enable the full shape optimization of plate structures, involving all geometrical variables. This is achieved by deriving the exact/implicit derivatives of the BEM plate formulations with respect to plate thickness for the first time.
Previous work by the research community involving exact⧵implicit derivatives of BEM formulations have focused on 2D structures, 1,3,21-28 with some work conducted on 3D structures, 13,29 and one work so far on plate structures. 20 One notable example is Huang et al. in which the implicit derivatives for the 2D dual boundary element method (DBEM), a version of BEM effective at modeling cracks, were developed for the first time and used to estimate the reliability of 2D structures using the First-Order Reliability Method (FORM). Another notable example is Brancati et al., 13 in which the implicit derivatives of 3D boundary element formulations were used to optimize noise levels in an aircraft cabin.
In summary, the main novelty of this work is that the exact/implicit derivatives of the BEM plate formulations with respect to plate thickness have been derived for the first time, enabling the full shape optimization of plate structures with the BEM. These implicit derivatives will be validated against derivatives obtained from the finite difference method (FDM) and from an analytical solution. To demonstrate the full shape optimization of plates with the BEM, a numerical example involving RBDO and manufacturing cost optimization is presented.
The methodology behind manufacturing cost estimation is presented in Section 2, followed by the methodology used for RBDO in Section 3. The implicit derivatives of the BEM plate formulations with respect to plate thickness are presented in Section 4, alongside the validation of these implicit derivatives. The numerical example involving RBDO and manufacturing cost estimation is shown in Section 5. The implicit derivatives of the BEM plate fundamental solutions with respect to plate thickness are presented in detail in the Appendix.

PARAMETRIC APPROACH TO MANUFACTURING COST ESTIMATION
The parametric cost estimation methodology used in this work is based on the parametric cost estimation methodology described in detail in appendix C of the NASA Cost Estimating Handbook. 30 A parametric model can be created to estimate the manufacturing cost of a structural component. The parametric model used in this work is a linear regression of the following form: where Cost is the manufacturing cost of a structural component, and [x 1 , x 2 , … x k ] are independent variables that influence the manufacturing cost, also known as cost drivers. [ 1 , 2 , … k ] are regression coefficients. The cost drivers could be the dimensions of a structural component, such as its length and width. The creation of this parametric model required the use of a historical database containing the details of similar structural components, and could be created by collecting the details of structural components from a range of similar aircraft models. In this work, this database is an artificial database created by the authors.
To minimize the effects of multicollinearity, a stepwise regression approach is taken when creating the regression model seen in Equation (1). The most important cost drivers are identified based on their Pearson correlations coefficient with respect to the manufacturing cost. The most important cost driver is used as the sole independent variable of the first iteration of stepwise regression procedure. In the following iterations, the next most important cost drivers are added to the model. If the addition of a driver noticeably improved the model, it is kept in the model, otherwise it is removed. This is repeated until all of the important cost drivers have been tried in the model. Therefore, the resulting regression model will include only the very most important cost drivers, mitigating the effects of multicollinearity and providing an accurate parametric model for estimating the manufacturing cost of a structure.

RELIABILITY-BASED DESIGN OPTIMIZATION (RBDO)
In the field of reliability analysis, the boundary between succeeding or failing to meet a certain set of criteria can be represented mathematically by a limit state function (LSF) g(Z). For example, if we are looking at the probability of a structure failing due to load, the LSF will be: where Z is a vector of random variables (Z ∈ R n r where n r is the number of random variables), and Q is a subset of Z if R is a random variable, where R is the resistance of the structure to some load S. If S(Q) > R then g(Z) < 0 and the structure is considered to have failed, while if S(Q) ≤ R then g(Z) ≥ 0 and the structure is considered safe. The probability that the set of criteria has failed to be met is termed the probability of failure P F , while the probability that the set of criteria has been successfully met is termed reliability P R . In the example outlined above, these probabilities would correspond to the probabilities of the structure failing or being safe under the load S respectively. Reliability can be determined by evaluating the following integral: where f Z (Z) is the joint PDF of Z. P R , and P F are obtained by integrating over the failure region (g(Z) < 0) and the safe region (g(Z) ≥ 0) respectively. All of the design variables are assumed to be mutually independent. The integral seen in Equation (3) can be difficult to evaluate if there are many variables in X or if the boundary g(Z) = 0 is non-linear. Therefore, several methods have been developed to evaluate the integral in Equation (3). The most widely known are Monte Carlo Simulations (MCS), the FORM, and the Second-Order Reliability Method. This work will focus on the FORM due to its efficiency. The reliability P R shown in (3) can be represented in terms of a reliability index as: where Φ denotes the CDF of the standard normal distribution. A large value for the reliability P R corresponds to a large value for the reliability index . can be found by rearranging the above equation to yield: where Φ −1 is the inverse of the CDF of the standard normal distribution. RBDO involves optimizing the design of a structure such that the reliability of the structure achieves a certain level of reliability. There are two main approaches to RBDO, the reliability index approach (RIA) and the performance measure approach (PMA).

Reliability index approach (RIA) to RBDO
In the RIA, the optimization problem is: where d is a vector of n d design variables (d ⊆ (Q) and n d ≤ n r ), FORM is the reliability index from FORM, and target is the target reliability index. d L and d U are vectors containing the lower and upper bounds respectively of the design variables. For a given vector of random variables Z, the RIA calculates the reliability index via FORM for each iteration of the optimization procedure.

Performance measure approach (PMA) to RBDO
In the PMA, the optimization problem is: where g(Z * ) is the value of the limit state function evaluated at the most probable point (MPP) found from the PMA. The PMA algorithm used in this work is the hybrid mean value (HMV) algorithm, 31 due to its enhanced efficiency and stability. The PMA can be thought of as the inverse of the RIA. For a given vector of random variables Z, the PMA calculates the MPP point Z * for which FORM = target . This MPP Z * is then used by the next iteration of the optimization procedure.
The PMA requires the derivatives of the constraints, that is, the limit state function g in Equation (7), to be calculated. In this work, g will be a function of boundary stresses or internal displacements. Therefore, the derivatives of the boundary stresses and internal displacements need to be derived. Since the BEM is used in this work to calculate boundary stresses or internal displacements, the derivatives of the BEM formulations for plate structures will need to be derived. These derivatives have been derived for the first time and are presented in the next section.

RBDO CONSTRAINT DERIVATIVES
To improve the computational efficiency of conducting RBDO for plate structures using the PMA approach, the exact derivatives of the response of a BEM plate model were derived for the first time.
In this work, Latin letter indexes (e.g., i, j, k) can take values from 1 to 3, while Greek letter indexes (e.g., , , , ) can take values of either 1 or 2.

BEM formulations for plates
In this section, the BEM formulations for plates are presented.

F I G U R E 2 Sign convection for displacement, rotations, tractions, and moments for shear deformable plates 33
Consider the plate of thickness h shown in Figure 1. The x 1 − x 2 plane is the middle surface x 3 = 0, or membrane, of the plate where −h∕2 ≤ x 3 ≤ +h∕2. x 3 can be described in terms of a non-dimensional variable x 3 , such that The displacements of the plate are u , where u 1 and u 2 are the displacements in the directions x 1 and x 2 respectively. The rotations of the plate are w , where w 1 and w 2 are the rotations of the plate in the directions x 1 and x 2 respectively. The displacement of the plate in the direction x 3 is w 3 . The tractions are denoted as t and p i . t are tractions due to membrane stress resultants N . p are tractions due to bending stress resultants M . p 3 denotes traction due to shear stress resultants Q . A diagram explaining the sign convention of these displacements, rotations, and tractions for plate bending can be seen in Figure 2

Boundary integral equations
From Reference 34, the discretized BEM boundary integral equations for plate bending are shown below. The discretized BEM boundary integral equation for the plate membrane is: where: where x n c , n c = 1, 2, … , N c (where N c = N n ) is the collocation node, N n is the number of nodes, N e is the number of elements, and M is the number of nodes per element. M = 3 in the case of quadratic elements. Superscript m indicates that an equation is for the plate membrane. S n e is the shape function of node of element n e , and J n e is the Jacobian of element n e . P m,n e and Q m,n e are fundamental solutions evaluated at node of element n e . u n e are in-plane displacements at node of element n e , and t n e are tractions due to membrane stress resultants N at node of element n e .
The discretized form of the displacement boundary integral equation for bending is: where: where: where superscript b indicates that an equation is for plate bending. denotes traction due to shear stress resultants Q , at node of element n e . q 3 is constant uniform loading over the entire top surface of the plate. is the Poisson's ratio, and = √ 10∕h is the shear factor, where h is the thickness of the plate. In the above equations, T b ij , U b ij , and V b i, are the fundamental solutions for plate bending, while T and U are the fundamental solutions for the membrane. Expressions for these fundamental solutions can be found in Reference 32. The integral symbol − ∫ represents Cauchy principal value integrals. w 1 and w 2 denote rotations in the directions x 1 and x 2 respectively, and w 3 denotes displacement in the direction x 3 . u 1 and u 2 are the displacements in the directions x 1 and x 2 respectively. p k are the bending and shear tractions with p = M n and p 3 = Q n . t 1 and t 2 are membrane tractions in the directions x 1 and x 2 respectively where t = N n . The integrations are carried out over the boundary S of the structure's domain. The terms C m and C b ij are free terms and their values can be directly evaluated from a consideration of rigid body motion. 32 Since the fundamental solutions shown in the integral Equations (9), (10), (49)-(52) are of the order of ln(1∕r) or 1∕r, (where r is the distance between the collocation node and the field point) mathematical singularities can occur when the collocation node lies within the same element as the field point. Weakly singular integrals are defined as integrals with singularities of the order ln(1∕r) or 1∕r such as those seen in Equations (49)-(52), (9) and (10). In this case, the transformation of variable technique proposed by Telles 35 is used. For Equations (49) and (9), rigid body motion is also applied. For each of the integral equations seen above, when the collocation node is near to the field point, but is not in the same element as the field point, the integral shows near-singular behavior. In this case, the element subdivision technique is used. Details on these methods can be found in Reference 32.
The system of equations used in the BEM is of the form Hu = Gt. Where H is a (5N n × 5N n ) matrix, and G is a (5N n × 5N e M) matrix. u is a (5N n × 1) vector of known and unknown displacements, and t is a (5N e M × 1) vector of known and unknown tractions. The final system of equations can be written as: where A is a (5N n × 5N n ) usually unsymmetric and dense coefficient matrix composed of parts of H and G, X is a (5N n × 1) vector containing all of the unknown boundary displacements and tractions, and F is a (5N n × 1) vector containing parts of G multiplied by known tractions and parts of H multiplied by known displacements. This system of equations, after applying the boundary conditions, can be solved using LU decomposition. An example of a BEM mesh for a simple plate structure can be seen in Figure 3. As seen in Figure 3, the boundary of a structure is discretized using continuous quadratic elements except at the corners-where due to the non-uniqueness of the normals, semi-discontinuous quadratic elements are used.

Displacements and rotations at internal points
The displacements and rotations at internal points can be evaluated using the solution X to the system of equations. The in-plane displacements at some internal point X n i , n i = 1, 2, … , N i (where N i is the number of internal points) are: where: The rotations w 1 and w 2 , and vertical displacement (out-of-plane displacement) w 3 , are: where: and:

Boundary stress resultants
In shear deformable plate theory, the normal stress components are assumed to vary linearly through the thickness of the plate. Therefore: where N are the membrane stress resultants, M are the bending stress resultants, and −h∕2 ≤ x 3 ≤ +h∕2 (see Figure 1). In this work, an indirect approach is used to evaluate boundary stresses. The boundary stresses are evaluated from boundary tractions and tangential strains. More detail on this method can be found in Reference 32. A brief description of this method is outlined below. A local coordinate system can be defined on a boundary element n e such that ê where ê n e = ê n e ( ), and x n e are the global coordinates of node of element n e . Therefore, the rotation matrix ê for node of element n e can be written as: This local coordinate system can be seen in Figure 4. .
The local boundary stresses can now be calculated using Equations (28) and (29). The global stresses N n e can be calculated from the local stressesN n e via: Using the relationships in Equation (27), the global membrane stress resultants can be obtained: The tangential component of the rotations are: The local tangential strain is:̂n .
The local moments and shear stresses are:M    .
The global moment stress resultants for node of element n e are:

Derivatives of BEM formulations for plates with respect to plate thickness h
In this section, the BEM formulations for plates with respect to plate thickness h are presented.

Boundary integral equations
The derivatives of the Boundary integral equations seen in Section 4.1 with respect to plate thickness h are presented here. The derivative of the discretized BEM boundary integral equation for the plate membrane (Equation 8) with respect to plate thickness h is:  ) with respect to plate thickness h is: where: where: are the derivatives of the free terms seen in Equations (8) and (11) respectively and their values can be directly evaluated from a consideration of rigid body motion.
In BEM-based IDM the system of equations is H ,h u + Hu ,h = G ,h t + Gt ,h , where H, G, u, and t are the same as defined in Section 4.1, and H ,h , G ,h , u ,h , and t ,h are their derivatives. This system of equations can be rewritten as: where A and X can be obtained from Equation (16). Since the right-hand side of Equation (53) is known, LU decomposition can be used to obtain the unknown derivatives of boundary displacements and tractions X ,h .

Displacements and rotations at internal points
The derivatives of the displacements and rotations at internal points seen in Section 4.1.2 with respect to plate thickness h are presented here. The derivatives of the in-plane displacements at some internal point X n i , n i = 1, 2, … , N i (where N i is the number of internal points) are: where: The derivatives of the rotations w 1 and w 2 , and vertical displacement (out-of-plane displacement) w 3 , with respect to plate thickness h are: where: and:

Boundary stress resultants
The derivative of the through thickness stress in a plate (Equation 25) with respect to plate thickness is: where N ,h is the derivative of the membrane stress resultants, M ,h is the derivative of the bending stress resultants, and x 3,h is the derivative of x 3 where −h∕2 ≤ x 3 ≤ +h∕2 (see Figure 1). As mentioned in Section 4.1, x 3 can be described in terms of a non-dimensional variable x 3 , such that The derivative of the normal component of the local stress is: The derivative of the tangential component of the local stress is: The derivatives of the out-of-plane stress resultants can be calculated in a similar manner as the in-plane stress resultants.
The derivative of the normal components of the tractions at node of element n e are: .
The global moment stress resultants for node of element n e are:

Validation
To validate the derived BEM formulations, results were compared with an analytical solution for a square simply supported plate subjected to uniform constant pressure presented in Reference 36. The square plate has edges of length a and thickness h, it is subjected to a uniform pressure q 3 . All four edges of the plate are simply supported such that w 3 = 0 along the edges. The plate can be seen in Figure 5, and details of the plate are shown in Table 1.
The maximum vertical deflection w max 3 will occur in the center of the plate. The solution given in Reference 36 for the square plate is: where D is the flexural stiffness: . The derivative of Equation (80) with respect to the thickness h is: where D ,h = 3D∕h. Equation (82) is the exact solution for w max,h . The derivative w max 3,h was also calculated via IDM and the FDM for a range of stepsizes, and at different values for h. The percentage differences between the exact solution for w max 3,h and the values of w max 3,h obtained from the IDM and the FDM can be seen in Table 2. For both the FDM and the IDM, a BEM model was created of the square plate consisting of 32 quadratic elements. It can be seen that the stepsize used with the FDM has a significant impact on the difference between the FDM and the exact solution for w max 3,h . Higher stepsizes cause greater instability, while lower stepsizes provide greater accuracy. The IDM results are in excellent agreement with the analytical solution.
In conclusion, the IDM has been successfully validated against an exact solution, and against the FDM. The IDM was shown to be more accurate than the FDM for most of the stepsizes tested. Furthermore, unlike the FDM, its accuracy was not dependent on the value of stepsizes, indicating the IDM is more robust and stable than the FDM.

NUMERICAL EXAMPLE
Section 4.3 demonstrated the high accuracy and robustness of the newly developed IDM. To now demonstrate its efficiency, the IDM is employed as part of a numerical example featuring the plate structure seen in Figure 6. In this numerical example, the geometrical design of this plate will be optimized in terms of its manufacturing cost and its reliability. The optimization procedures used in this numerical example require the calculation of constraint derivatives, which can be calculated using either the IDM or the FDM. To determine the efficiency of the IDM, a comparison is made between the time required by IDM and the FDM to complete the optimization at the end of this numerical example. The plate is simply supported around its outer edge such that the vertical deflection of the plate is zero on its outer edge, and it is subjected to boundary tension and bending moments along this outer edge. It is composed of Aluminum 6061-T6, an aluminium alloy commonly found in aircraft structures due to its high strength and low weight. It has a Young's modulus E of 68.9 GPa, a Poisson's ratio of 0.33, and a tensile yield strength of 276 MPa. Material properties of Aluminum 6061-T6 can be found in Reference 37. The geometry of the plate is described by the geometrical variables: W 1 , L 1 , R 1 , W 2 , L 2 , R 2 , and h. These geometrical variables, along with the boundary tractions, boundary bending moments, and yield strength, are treated as random variables in the optimization procedure. Details of the random variables can be seen in Table 3. A total of 112 quadratic elements were used in the BEM mesh of this plate structure.
The manufacturing cost of the plate is estimated using the approach outlined in Section 2. To use this approach, a database was created that contains the details of 100 plates with a geometry similar to the plate shown in Figure 6. Details  of this database can be seen in Table 4. By using this database with the approach outlined in Section 2, the following formula for the manufacturing cost of the plate can be obtained: When run with the database, this formula demonstrated a high coefficient of determination of R 2 = 0.82 and R adj = 0.81, indicating the high level of accuracy associated with this formula. By using this formula, the optimization of the plate structure can now be conducted with respect to manufacturing cost.
To investigate the impact of optimizing the plate with respect to reliability, two different approaches are investigated for the optimization of the plate structure: • Approach 1: Both the manufacturing cost and the reliability of the plate are considered during the optimization procedure.
• Approach 2: Only the manufacturing cost of the plate is considered during the optimization procedure.
TA B L E 4 Details of the database containing 100 plates with a geometry similar to the plate shown in Figure 6 Variable Both approaches investigated in the following sections. It is expected that the optimal designs obtained from Approaches 1 and 2 will be significantly different.

Manufacturing cost optimization and RBDO
In this section, both the manufacturing cost and the reliability of the plate are considered during the optimization procedure. In this case, the optimization problem is: where Z = [W 1 , L 1 , R 1 , W 2 , L 2 , R 2 , h, t, p, y ] is a vector of the random variables from Table 3 that influence g, and Q = [ is a vector of the random variables from Table 3 that influence the maximum Von-Mises stress max in the plate structure (Q is identical to Z, except Q lacks y ). The IDM, with the PMA outlined in Section 3.2, is used to solve the optimization problem seen in Equation (84). Similar to the previous section, the optimization was conducted using Matlab © with the nonlinear multivariable optimization routine "fmincon." The optimal plate designs from "fmincon" for a range of target reliability indices target can be seen in Table 5, and diagrams of these optimal designs can be seen in Figure 7. It can be seen from Table 5 and Figure 7 that the plate becomes thicker and wider as target increases, which makes sense given that the maximum stress is expected to occur at the corners of the central hole. To verify that these optimal designs achieve the desired reliabilities, the reliability index from the FORM FORM was calculated for each of the optimal designs and are presented in Table 5. It can be seen that FORM and target are in very good agreement; the difference between them is less than 0.1%.
The target reliability indices target and the corresponding target probabilities of failure P F,target for each of the optimal designs can be seen in Figure 8. It can be seen that the probability of failure drops by a magnitude of 8 as the reliability index increases; the probability of failure corresponding to a reliability index of 2 is 2.3%, while for a reliability index of 6 it is 9.9 × 10 −8 %. The cost on the other hand, increases from 24.92 € for a reliability index of 2, to 81.75 € for a reliability index of 6, a percentage increase of 236%. The optimization history from the IDM for a range of target reliability indices can be seen in Figure 9. It can be seen that the significantly more iterations were required to obtain convergence for the case with target = 6. This can partly be explained by the fact that the jump in the probability of failure between target = 5 to target = 6 was by a magnitude of 3, while previous jumps in the probability of failure were only by a magnitude of 2 or 1. Therefore, significantly more alterations were required to the geometry of the design, increasing the number of iterations needed.
Up until this point, only the optimization results with the IDM have been presented. This is because Section 4.3 demonstrated the high accuracy and robustness of the IDM, and its excellent agreement with the FDM when calculating the constraint derivatives required by the optimization procedure. To demonstrate the efficiency of the IDM, the CPU time required by the IDM to optimize the plate structure was compared to the CPU time required by the FDM. The average CPU time required by the IDM and the FDM to complete one optimization iteration was calculated by averaging over 100 optimization iterations. The results are shown in Table 6. It was found that the IDM was, on average, 19% faster than the FDM. For a single iteration or for a simple structure, this increase in efficiency may not be significant, but it could be very useful for a more complex structure that requires many iterations. In this case, the IDM could significantly reduce the total optimization time.
A flowchart for designing the optimization code used in this section can be seen in Figure 10.

Manufacturing cost optimization
In this section, only the manufacturing cost of the plate is considered during the optimization procedure. In this case, the optimization problem is: Minimize where d = [W 1 , L 1 , R 1 , W 2 , L 2 , R 2 , h] is the vector of design variables, and n d = 7 is the number of design variables. The initial design of the plate d 0 , as well as the lower and upper bounds of the design variables, d L and d U , are the same as in the previous section. The optimization problem in Equation (86) was conducted using Matlab © with the nonlinear multivariable optimization routine "fmincon." The optimal plate design from "fmincon" can be seen in Table 7, and a diagram of this optimal design can be seen in Figure 11. It can be seen that the values of the design variables are equal to their lower or upper bounds, except for R 2 which stayed at its initial value since it is not included in the regression model seen in Equation (83). The manufacturing cost was reduced from 61.40 € for the initial design, to only 15.06 € for the optimal design; this represents a significant reduction in cost. However, when evaluating the reliability of this design using the limit state function shown in Equation (85), a probability of failure of P F = 99.6% is obtained. This probability of failure is significantly higher than that the highest obtained from the previous section, which was only 2.3%. Such a high probability is unacceptable, and highlights the importance of taking into account both manufacturing costs and reliability when optimizing the design of a structure.

F I G U R E 11
The optimal plate design when only the manufacturing cost of the plate is considered during the optimization procedure F I G U R E 12 A flowchart for designing the optimization code used in Section 5.2 A flowchart for designing the optimization code used in this section can be seen in Figure 12.

CONCLUSIONS
In conclusion, this article presented a novel methodology for the manufacturing cost and reliability-based optimization of plate structures with the BEM, with the aim of improving the accuracy, robustness, and efficiency of the optimization of plate structures. The derivatives of the BEM plate formulations, with respect to plate thickness, were derived for the first time and used as part of an IDM, enabling the full shape optimization of plate structures with the BEM. These implicit derivatives were validated against derivatives obtained from the FDM and from an analytical solution. The IDM was found to be in excellent agreement with the analytical solution, and more robust and accurate than the FDM for most of the step-sizes investigated. To demonstrate the efficiency of the newly developed IDM, it was employed as part of a numerical example involving the RBDO and manufacturing cost optimization of a plate structure. The design parameters in the optimization included all the geometric parameters describing the shape of the structure. Results indicate that the IDM is 19% faster, on average, in terms of CPU time than the FDM when performing this optimization. Such results represent a significant reduction in the computation time associated with the optimization of complex structures.

ACKNOWLEDGMENTS
This project has received funding from the Clean Sky 2 Joint Undertaking (JU) under grant agreement No. 864154 Project MASCOT. The JU receives support from the European Union's Horizon 2020 research and innovation programme and the Clean Sky 2 JU members other than the Union.

CONFLICT OF INTEREST
The authors declare that they have no conflict of interest.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.