Gradient-based algorithms are deterministic methods that can significantly reduce the computation time required to converge to an optimal solution by limiting search to promising directions.
Gradient based approach algorithms have been compared to genetic algorithms in terms of calculation time and the choice of objective functions. They are mainly used because of their speed and however they are very sensitive to the initial condition [81] and in this sense they are not robust [341]. On the other hand, gradient based algorithms can lack in global optimality but they allow multiple constraints, which can be very useful for complex engineering designs. More often for complicated problems, it is difficult to obtain a global optimal because conventional algorithms (such as feasible direction methods) are susceptible to converge to the local optimal point [342]. Therefore the user is prompted to interfere in the design process and adjust the design parameters or shift the initial feasible domain. For example, Fuglsang et al. [34] apply a Sequential Linear Programming (SLP) [343] when the design vector was feasible and the Method of Feasible Directions (MFD) [344] was used to return the design vector to the feasible domain when it was unfeasible. Kenway et Martins [53] use SNOPT [345], an optimizer based on the Sequential Quadratic Programing (SQP) approach. Similarly, Ning et al. [36] use central differencing and a multi-start approach to improve the convergence behavior of the SQP algorithm.
Hybrid methods between GA and GBA have been investigated by [81,88,319,346,347]. Grasso [88] implement this scheme by first using the GA to explore large domain that contain less local optima problems, and an optimal solution is found. This latter is used as the initial configuration for the GBA which searches for an accurate optimal solution in the smaller design space. Bizzarrini [81] compares the results of a hybrid scheme and a single-algorithm (GA) and shows that the hybrid is more effective with higher accuracy and low sensitivity to local minima.
3 Application of gradient-based methods to aerodynamic optimisations
Gradient-based optimisation is a calculus-based point-by-point technique that relies on the gradient (derivative) information of the objective function with respect to a number of independent variables. The nature in which gradient-based methods (GBM) operate make them well suited to finding locally optimal solutions but may struggle to find the global optimal [68]. With gradient-based algorithms an understanding of the design space is assumed, as an appropriately pre-conceived starting design point must be given. Kenway and Martins [13] point out that with increasingly higher fidelity aerodynamic optimisations, a more refined initial design should be used so that the optimisation does not diverge too far from the baseline. If large changes in topology are expected lower fidelity panel codes, such as that described by Vecchia and Nicolosi [27], can facilitate useful optimisation procedures. Typically, the higher the fidelity analysis used the more compact the design variables will need to be to allow effective optimisation with a gradient based optimiser.
Gradient-based optimisation is, in its most basic form, a two step iterative process which can be summarised mathematically as:
(1)
where ∇f is the gradient of function F(X), and X is a vector of the design variables. The first step is to identify a search direction (gradient), ∇f, in which to move. The second step is to perform a one-dimensional line search to determine a distance/step size h along ∇f that achieves an adequate reduction of some cost function, i.e. define how far to move in the search direction until no more progress can be made [69]. A schematic diagram illustrating the operation of a gradient-based optimisation is shown in Fig. 9. In-depth benchmarking of gradient based algorithms for aerodynamic problems has been conducted by Secanell and Suleman [70] and Lyu et al. [71].
Fig. 9. Schematic diagram of a gradient-based aerodynamic optimisation process.
Gradient based algorithms are extensively used in aerospace optimisation as they exhibit low computational demands when handing many hundreds of design variables – this makes them well suited for optimising shapes based on deformative geometric parametrisations [12,64]. Significant difficulties arise if they are not applied within a restricted set of functions with well defined slope values due to a dependency upon the existence of derivative information via some sensitivity analysis. There are several different methods for sensitivity analysis for which four general classes can be distinguished: 1) finite-difference methods; 2) complex-step derivative approximation; 3) automatic/algorithmic differentiation; and 4) analytic methods. It is important to understand their relative merits since none are a clear choice for all classes of problem. Comparative studies on the numerical sensitivity analysis for aerodynamic optimisation has been conducted by Martins et al. [72] and Peter and Dwight [73], while Martins and Hwang [74] offer a detailed discussion for computing derivatives within multi-disciplinary computational models.
The computational expense of evaluating gradients using finite-difference or the complex-step method provide a simple and flexible means of estimating gradient information, but are considered excessive with respect to hundreds of variables [75]. These approaches preserve discipline feasibility, but they are costly and can be unreliable. Finite-differencing, while not used to provide gradients for the optimisation itself have been used by Kenway and Martins [13] and Skinner and Zare-Behtash [67] to provide gradients for stability derivative constraints. For a restricted number of design variables the complex-step method is suitable for sensitivity analysis as demonstrated by Kenway and Martins [76]. They employ the complex step method to provided gradients for the Sparse Nonlinear Optimiser (SNOPT) algorithm, originally developed by Gill et al. [77], in the constrained optimisation of wind turbine blades. It is commented on that to increase the dimensionality of the problem an analytic sensitivity analysis would have to be adopted.
Finite-differencing or complex-step methods employed for providing sensitivity analysis for low-fidelity codes can be considered appropriate due to low computational demand. Ning and Kroo [66] optimise a series of wing topologies investigating fundamental wing design trade-offs for which sensitivity analysis of the objective and constraints were approximated by finite-differencing. Results provided by the sequential quadratic programming method show robust and quick convergence able to determine relative gradients between approximated area-dependant weight, effects of critical structural loading, and stall speed constraints.
In the presence of several hundred design variables and constraints the analysis code will require a particularly long time to evaluate sensitivities. Automatic/algorithmic differentiation or analytic derivative calculations (direct or adjoint) can be used to avoid multi-discipline analysis evaluations. Pironneau [78], pioneered the adjoint method in fluid dynamics, showing that the cost of computing sensitivity information was almost completely independent of the number of design variables, and hence the overall cost of optimisation is roughly linearly proportional to the number of design variables. Lyu et al. [71] more recently demonstrated that both SNOPT and sequential least squares programming [79,80] (SLSQP) gradient-based algorithms with analytically derived adjoint gradients require far fewer total functions calls when compared to using finite-differencing for high-fidelity large scale aerodynamic optimisations. The adjoint form of the sensitivity information is particularly efficient for aerodynamic optimisation applications as the number of cost functions (outputs) is small, while the number of design variables (inputs) is relatively larger.
The discrete adjoint method (as opposed to continuous adjoint method) is generally favoured in aerospace-based optimisation as it ensures that sensitivities are exact with respect to the discretised objective function [81,82]. The implementation of the adjoint method for the governing equations of the flow analysis can often be difficult to derive and require direct manipulation; adjoint methods require much more involved detailed knowledge of the computational domain. One way to approach this difficulty is to use automatic/algorithmic differentiation, which is a method based on the systematic application of the differentiation chain rule to the source code to compute the partial derivatives required by the adjoint method. Mader et al. [83] developed a discrete adjoint method for Euler equations using automatic differentiation, later followed by Lyu et al. [84] who extended and developed this adjoint implementation to Reynolds-averaged Navier-Stokes (RANS) equations and introduced simplifications to the automatic differentiation approach. Methods developed have shown robust an efficient application to high-fidelity optimisation [12,63].
Hicken and Zingg [85] adopted similar methods for the high-fidelity aerodynamic optimisation of non-planar wings addressing the non-linearity of wake shape and how it can impact the induced drag. Several non-planar geometries, inherently creating non-planar wake-wing interactions, are optimised using discrete adjoint sensitivities and the SNOPT algorithm. This work illustrates the drawbacks in static-wake assumptions, demonstrating that higher-order effects must be included for accurate induced drag prediction and hence for meaningful optimisations. This work was followed by Gagnon and Hicken [86] for the aerodynamic optimisation of un-conventional aircraft configurations; adapted optimisation results are shown in Fig. 10. Here, for the aerodynamic metrics, the gradients are evaluated analytically by using the discrete-adjoint variables while other gradients are provided by the complex-step method. This work is notable as it enables used axial deformation combined with free-form deformation to achieve both local and global geometric manipulations, the effects of this can be seen in Fig. 10, allowing span, sweep, dihedral, taper, twist and aerofoil sectional shape changes. These global geometric variables are generally not considered in high-fidelity simulation. The cost of allowing such geometric variation away from the baseline under high-fidelity optimisation limited how many variables could be considered in any one optimisation process. The authors observed limited optimisation in some wing configurations because of this.
Fig. 10. Upper surface pressure coefficient contours over initial and optimised un-conventional aircraft configurations [86].
Aeroelastic optimisation requires the coupling of aerodynamic and structural models for most effective sensitivity analysis in optimisation routines. Even small changes in aerodynamic shape can have a large influence on aerodynamic performance with various flow conditions resulting in multiple shapes. Wing flexibility impacts not only the static flying shape but also it's dynamics, resulting in aeroelastic phenomenon such as flutter and aileron reversal. Based on this principle, to enable high-fidelity aerostructural optimisation while encompassing hundreds of design variables, Martins et al. [87] proposed the use of a coupled adjoint method to compute sensitivities with respect to both the aerodynamic shape and the structural sizing. Kenway et al. [13,88] subsequently made several developments and demonstrated that the computation of coupled aeroelastic gradient calculations were scalable to thousands of design variables and millions of degrees of freedom, and since applied it to the aerostructural optimisation of high aspect ratio wings with different structural properties [89]. More recently, Burdette et al. [90,64] applied the coupled discrete adjoint method with the sparse non-linear optimiser SNOPT for wing morphology optimisation. This approach was capable of handling over a thousand design variables and constraints. The coupled adjoint method is also applicable to lower-fidelity models in studies like that by Elahm and Tooren [91], where they used a vortex lattice method and finite element analysis tool capable of accurately mimicking high-fidelity accuracy at a greatly reduced computational cost.
The coupling of design constraints makes the optimiser additionally capable of considering more sophisticated criteria. Mader and Martins [92] included flight dynamics into the coupled adjoint sensitivity and explored the use of static and dynamic stability constraints. Their result showed that coupling stability constraint sensitivities into the adjoint formulation had a significant impact on optimal wing shape. Elsewhere, structural dynamics were considered by Zhang et al. [54] who used a coupled-adjoint formulation to include flutter constraints. The flutter constraints used the coupled aerodynamic/structural solver to suppress flutter onset by identifying dominant modes and adjusting variables such as the wing stiffness.
Grossman et al. [93] investigated using modular sensitivity analysis for aerostructural sequential optimisation of a sailplane. They showed that coupled aerostructural optimisation gave higher performance designs than those identified by sequential optimisation of aerodynamics followed by structural optimisation. Subsequently, Grossman et al. [94] optimised the performance of a subsonic wing configuration showing that while modular sensitivity analysis for sequential optimisation reduced the total number of function calls and sensitivity calculations, the wing performance gain was limited. When performing sequential optimisation the optimiser does not have sufficient information necessary for aeroelastic tailoring. This limitation of sequential optimisation is further explained by Chittick and Martins [95].
A significant drawback of all gradient-based algorithms is the requirement for continuity and low-modality throughout the design space otherwise the algorithm may become sub-optimally trapped. The challenge is that an aerodynamic shape analysis throughout a geometrically varying search space will encounter both non-continuous topological and local flow changes, each providing local optima [96,71,97]. Gradient-dependant algorithms’ robustness significantly decreases in the presence of discontinuity and lack of convergence, usually related to turbulence modelling, making the objective function noisy [98]. Kenway [99] encountered such a problem with aerodynamic shape optimisation with a separation-based constraint formulation to mitigate buffet-onset behaviour at a series of operating conditions. The discontinuity from the ‘separation sensor’ function arose from monitoring the wing local surface for separated flow; this resulting in locally negative skin friction coefficients. To address this issue blending functions were to be implemented to smooth the discontinuity, smearing the separation sensor value around the separated flow region.
Kenway and Martins [99], among several others [100,101,64,13], have used multi-point optimisation strategies in order to consider several operating conditions simultaneously. For more realistic and robust design it is crucial to take into account more than one operating condition, especially off-design conditions, which form additional multi-objective requirements into the optimisation. Fig. 11 shows the results for both single-point and multi-point optimisation of the Common Research Model configuration presented by Kenway and Martins; [99] results are shown for the nominal operating condition. The single-point optimisation achieved an 8.6 drag count reduction and the shock-wave over the upper surface of the wing is almost entirely eliminated. Drag divergence curves in this work show the nature of the single-point optimisation presenting a significant dip in the drag at the design condition, but the performance is significantly deteriorated at off-design conditions relative to the baseline condition. The multi-point optimisation, accounting for 3 design conditions, found that drag at the nominal operating condition increased by 2.8 counts and produced double shocks on the upper surface of the wing visible in Fig. 11. However, at the sacrifice of performance at the nominal operating condition, off-design conditions for the multi-point optimisation design was found to perform substantially better over the entire range of Mach numbers. Though biasing the optimisation toward certain operating conditions the authors show that multi-point optimisation with all conditions near the on-design condition is not sufficient for an overall robust design when considering operational envelopes.
Fig. 11. High performance low drag solutions found for single design point at nominal operating conditions. For multi-point optimisation performance at the nominal operating condition is sacrificed [99].
In order to determine the optimal device placement, a gradient-based optimization procedure based upon the optimality criteria and relevant performance sensitivities can be used (Takewaki et al., 2012). For a better understanding, various steps of this procedure are illustrated schematically in Fig. 5.51. In such a procedure, sensitivity analyses of the objective function should be usually performed with respect to the dampers’ damping coefficient to find the highest performance sensitivity. This method is subject to iterations (i.e., repeating the sensitivity analyses and finding the highest performance sensitivity) until the required total amount of supplemental damping is obtained. This category of algorithms usually requires programming.
Figure 5.51. Representative schematic diagram of gradient-based optimization procedures.
Adapted from Takewaki, I., 2009. Building Control with Passive Dampers: Optimal Performance-Based Design for Earthquakes. John Wiley & Sons (Asia), Singapore.
Typical gradient-based algorithms are often relative to linear structural behavior. Several studies have been conducted in the past for different device categories, such as:
•
Viscous dampers (e.g., Balling and Pister, 1983; Takewaki, 1997, 1999, 2000; Singh and Moreschi, 2001; Mahendra and Moreschi, 2001; Uetani et al., 2003; Lee et al., 2004; Lavan and Levy, 2005, 2006; Attard, 2007; Aydin et al., 2007; Cimellaro, 2007; Viola and Guidi, 2009; Aydin, 2012; Adachi et al., 2013; Lavan, 2015)
•
Viscoelastic dampers (Hwang et al., 1995; Singh and Moreschi, 2001; Mahendra and Moreschi, 2001; Lee et al., 2004; Park et al., 2004; Fujita et al., 2010)
•
Hysteretic dampers (e.g., Uetani et al., 2003)
•
TMDs (e.g., Zuo and Nyfeh, 2004; Wang et al., 2009; Salvi and Rizzi, 2014)
The search for an optimal state is one of the most fundamental principles in our world. Optimization is to find the best solution to a certain designated problem. Numerical methods for optimizing the performance of engineering problems have been studied for many years. For optimization of an objective, different categories of optimization techniques namely “Gradient based” methods and “Non-gradient based” methods can be used. The basic idea behind the gradient based methods is that a given function reaches its extremes (minimum or maximum) in the direction of its gradient. Gradient based methods are in general computationally faster (they require fewer objective function evaluations in case of problems with low number of design variables) than non-gradient based methods. Their main drawback is that they might converge to local minima, and their convergence to global minima depends mainly on the chosen starting point (initial guess) by the user (Fig. 7a). On the other hand, non-gradient based methods, do not rely, most of the time, on strong mathematical basis and make use of neither the gradient or the second derivative of the objective function as a direction of descent. These methods work based on function evaluations alone. In principle, they attempt to mimic nature in order to find the optimum of the objective function. One of the key features of these algorithms is that they search from multiple points in the design space, instead of moving from a single point like gradient based methods. Thus, although in general there is no proof that these methods converge to global optima, experience entails that they converge to global optima in most cases (Fig. 7b).
Fig. 7. Optimization methods: (a) Gradient based (b) Non-gradient based.
The main drawback of non-gradient based methods is that these algorithms are generally computationally slower than the gradient based ones. Indeed, the non-gradient based algorithms can require thousands of function evaluations and, in some cases, become non-practical. In order to overcome these difficulties, the so-called Hybrid algorithms, which take advantage of the robustness of the non-gradient based methods and the fast convergence of the gradient based methods, have been proposed by different scholars. A set of analytically formulated rules and switching criteria can be coded into the program to automatically switch back and forth among the different algorithms as the iterative process advances. Each technique provides a unique approach with varying degrees of convergence, reliability and robustness at different stages during the iterative optimization process.
It is also interesting to introduce to the reader the multi-objective optimization problems vs. single-objective optimization problems. In a large number of problems, there exists a need to find optimal solutions due to more than one objective. In this case, a multi-objective approach must be employed. Engineering problems demanding low cost, high performance and low losses are an example of applications where this approach is needed. For the single-objective optimization problems, a unique optimal solution exists. However, for multi-objective optimization problems, there exist a set of compromised solutions, known as the Pareto-optimal solutions or non-dominated solutions, which are based on the competing objectives. The goal of multi-objective optimization is to find such set of solutions. Once the solutions are obtained, the designer can choose a final design with further considerations. Non-gradient based methods are more suitable for multi-objective optimization problems since they are able to find the entire range of Pareto-optimal solutions. Gradient based methods typically use a weighting method to combine different objectives into one for handling multi-objective optimization problems. In this case, results would be dependent on the chosen weighting coefficients.
An overview of different gradient based and non-gradient based optimization algorithms are provided in Appendix A. Clearly, knowledge about the nature of the problem is a requirement for choosing the most suitable optimization tool for an application.
An extensive number of optimization methods based on algorithms have been presented in the literature to identify optimal parameter and placement of dynamic modification devices. In general, the most common ones can be categorized as gradient-based algorithms; GAs; control theory–based algorithms; and heuristic algorithms. The major applications of these optimization algorithms have been related to distributed-type damping systems, especially viscous dampers. In the following, the different procedures are briefly reviewed. Interested readers should refer to the references provided.
5.6.2.1 Gradient-based algorithms
In order to determine the optimal device placement, a gradient-based optimization procedure based upon the optimality criteria and relevant performance sensitivities can be used (Takewaki et al., 2012). For a better understanding, various steps of this procedure are illustrated schematically in Fig. 5.51. In such a procedure, sensitivity analyses of the objective function should be usually performed with respect to the dampers’ damping coefficient to find the highest performance sensitivity. This method is subject to iterations (i.e., repeating the sensitivity analyses and finding the highest performance sensitivity) until the required total amount of supplemental damping is obtained. This category of algorithms usually requires programming.
Figure 5.51. Representative schematic diagram of gradient-based optimization procedures.
Adapted from Takewaki, I., 2009. Building Control with Passive Dampers: Optimal Performance-Based Design for Earthquakes. John Wiley & Sons (Asia), Singapore.
Typical gradient-based algorithms are often relative to linear structural behavior. Several studies have been conducted in the past for different device categories, such as:
•
Viscous dampers (e.g., Balling and Pister, 1983; Takewaki, 1997, 1999, 2000; Singh and Moreschi, 2001; Mahendra and Moreschi, 2001; Uetani et al., 2003; Lee et al., 2004; Lavan and Levy, 2005, 2006; Attard, 2007; Aydin et al., 2007; Cimellaro, 2007; Viola and Guidi, 2009; Aydin, 2012; Adachi et al., 2013; Lavan, 2015)
•
Viscoelastic dampers (Hwang et al., 1995; Singh and Moreschi, 2001; Mahendra and Moreschi, 2001; Lee et al., 2004; Park et al., 2004; Fujita et al., 2010)
•
Hysteretic dampers (e.g., Uetani et al., 2003)
•
TMDs (e.g., Zuo and Nyfeh, 2004; Wang et al., 2009; Salvi and Rizzi, 2014)
•
TLCDs (Taflanidis et al., 2007)
5.6.2.2 Genetic algorithms
GAs are usually efficient to identify the optimal location of dynamic modification devices in buildings and they do not restrict to have a linear behaving structure. Such methods are found more suitable for problems where the PI is not a continuous function of the design variables (e.g., damping constant) (Singh and Moreschi, 2002). One of the advantages of this approach is the possibility to set a given dynamic modification device capacity, usually based on commercially available solutions. Consequently, the optimal location of a certain number of dynamic modification devices, with the fixed capacity selected, can be identified (Singh and Moreschi, 2002). As a shortcoming, these methods usually require long computational time demand (Singh and Moreschi, 2002).
The GA-based optimization techniques have been employed for the identification of optimal placement of various dynamic modification devices, such as:
•
Viscous dampers (e.g., Furuya et al., 1998; Singh and Moreschi, 2002; Wongprasert and Symans, 2004; Tan et al., 2005; Dargush and Sant, 2005; Silvestri and Trombetti, 2007; Lavan and Dargush, 2009; Kargahi and Ekwueme, 2009; Shin, 2010; Apostolakis and Dargush, 2010; Hejazi et al., 2013; Greco et al., 2016)
•
Viscoelastic dampers (e.g., Singh and Moreschi, 2002; Park and Koh, 2004; Dargush and Sant, 2005; Movaffaghi and Friberg, 2006; Shin, 2010; Qu and Li, 2012)
•
Metallic dampers (Dargush and Sant, 2005; Shin, 2010)
•
Hysteretic dampers (Moreschi and Singh, 2003; Ok et al., 2008)
•
Friction dampers (Moreschi and Singh, 2003; Miguel et al., 2014)
•
BRBs (Farhat et al., 2009)
•
TMDs (e.g., Hadi and Arfiadi, 1998; Arfiadi and Hadi, 2011; Singh et al., 2002; Marano et al., 2010; Mohebbi and Joghataie, 2011; Fu et al., 2011; Huo et al., 2013; Hervé Poh’sié et al., 2015; Venanzi, 2015; Greco et al., 2016)
•
Tuned liquid dampers (Ahadi et al., 2012; Chakraborty and Debbarma, 2016)
•
Isolation systems (Pourzeynali and Zarif, 2008; Charmpis et al., 2012; Xu et al., 2013)
5.6.2.3 Control theory–based algorithms
Control theory–based algorithms are advantageous due to their capability in reducing the computational efforts because they do not require to compute the structural response using dynamic analyses. The simplified sequential search algorithm (Lopez-Garcia, 2001; Garcia and Soong, 2002) and analysis-redesign procedure (Levy and Lavan, 2006) are among the simplest examples of this kind of algorithms. As noted by Whittle et al. (2012), a disadvantage of such methods may be the lack of PBD criteria within the methods. The control theory–based algorithms have been frequently used for viscous dampers (Zhang and Soong, 1992; Gluck et al., 1996; Lopez-Garcia, 2001; Yang et al., 2002a, bYang et al., 2002aYang et al., 2002b; Ribakov and Reinhorn, 2003; Main and Krenk, 2005; Levy and Lavan, 2006; Cimellaro and Retamales, 2007; Aguirre et al., 2013) and viscoelastic dampers (Zhang and Soong, 1992; Loh et al., 2000; Lopez-Garcia, 2001).
5.6.2.4 Heuristic algorithms
In addition to aforementioned algorithms, some other heuristic ones are applied in the literature for the purpose of optimizing dynamic modification systems, such as:
•
Viscous dampers: performance-based heuristic approach (Liu et al., 2005) and artificial bee colony algorithm combined with a gradient-based algorithm (Snomez et al., 2013)
•
Hysteretic dampers: adaptive smoothing algorithm (Murakami et al., 2013).
•
Friction dampers: backtracking search optimization algorithm (Miguel et al., 2015).
•
TMDs: bionic algorithm (Steinbuch, 2011), particle swarm optimization (Leung and Zhang, 2009), harmony search method (Bekdas and Nigdeli, 2011, 2013), ant colony optimization (Farshidianfar and Soheili, 2013), evolutionary operation (Islam and Ahsan, 2012), and charged system search (Kaveh et al., 2015)
As discussed in Chapter 3, numerical optimization techniques can be categorized as gradient-based and nongradient algorithms. Gradient-based algorithms often lead to a local optimum. Nongradient algorithms usually converge to a global optimum, but they require a substantial amount of function evaluations. For large-scale problems, as are often encountered in engineering design, nongradient algorithms are less desirable. Gradient-based algorithms require gradient or sensitivity information, in addition to function evaluations, to determine adequate search directions for better designs during optimization iterations.
In optimization problems, the objective and constraint functions are often called performance measures. Sensitivity, sensitivity coefficient, and gradient are terms used interchangeably to define the rate of change in a performance measure with respect to the change in design variable. They support not only gradient-based optimization but reveal critical design information that could guide designers to achieve improved designs in just one or two design iterations—for example, through the interactive design steps discussed in Chapter 3.
As shown in Figure 4.1, a change in the height h or width w of a cantilever beam with the rectangular cross-section affects the maximum bending stress inside the beam. The bending stress can be formulated as
FIGURE 4.1. Cantilever beam with a rectangular cross-section.
(4.1)
where P is a point force applied at the tip of the beam. The gradient of the bending stress can be calculated by taking derivatives of the stress in Eq. 4.1 with respect to height and width of the beam cross-section, respectively, as
(4.1a)
and
(4.1b)
From Eqs 4.1a and 4.1b, both gradients are negative, implying that increasing either height or width reduces the bending stress. Physically, it is obvious that increasing either height or width increases the moment of inertia I of the beam cross-section. In this case, I = wh3/12, therefore reducing the bending stress. Between the two design variables, increasing the height dimension h is two times more effective than the width w in reducing the bending stress if w = h (a square cross-section at the current design). Equations, such as Eqs 4.1a and 4.1b, provide desired information that supports designer to achieve improved design effectively.
The formulation of DSA can vary significantly depending on what type of design variables are being considered. In general, a structure consists of bars (or trusses), beams, membranes, shells, and/or elastic solid structural components. Depending on the constituents of the structure being designed, there are five different types of design variables—material, sizing, configuration, shape, and topology, as illustrated in Figure 4.2. For example, for a bar or beam structure shown in Figure 4.2a, material design variables can be the mass density ρ or modulus of elasticity E, while sizing design variables can be the cross-sectional areas A of individual bar members or the area moment of inertias I of the beam members. Configuration design variables are related to the orientations of components in built-up structures, such as the change shown in Figure 4.2b, in which orientation angles of individual members are altered in addition to length changes. Shape design variables describe the change in the length of one-dimensional (1-D) structures or the geometric shape of two-dimensional (2-D) and three-dimensional (3-D) structures, as illustrated in Figure 4.2c. Topology optimization determines the layout of the structure, as depicted in Figure 4.2d.
FIGURE 4.2. Illustration of design variables in different categories. (a) A built-up structure in which material parameters and cross-sectional areas of the bar members can be changed. (b) Configuration design by adjusting the orientations and lengths of truss members (Twu and Choi 1993). (c) Shape design for a 2-D engine connecting rod (Edke and Chang 2011). (d) Topology optimization of a solid beam (Tang and Chang 2001).
To carry out sensitivity analysis, it is obvious that the performance measure is presumed to be a differentiable function of the design, at least in the neighborhood of the current design. In general for structural designs, essential structural responses, such as displacement, stress, buckling load, frequency, and so forth, are differentiable functions with respect to design.
Sensitivity analysis for performance measures expressed explicitly in design variables, such as in Eq. 4.1, is straightforward. One derivative leads to sensitivity equations such as those in Eqs 4.1a and 4.1b. No expensive computation is required.
However, in most cases, performance measures cannot be expressed explicitly in design variables. In fact, in many cases, function evaluations must be carried out numerically using FEA. In these cases, there is no equation to take derivatives for. So, how do we calculate gradients to explore viable design alternatives and support batch mode design optimization?
In this chapter, we discuss sensitivity analysis for sizing and shape designs, as well as topology optimization. We use simple and analytical examples to illustrate the sensitivity analysis concept and detailed solution techniques. We start by introducing the basic formulation of a simple bar structure in Section 4.2, which serves as an example problem to facilitate our follow-up discussion. We then provide an overview on sensitivity analysis methods in Section 4.3. In Section 4.4, we discuss sensitivity analysis for sizing and material design variables, especially using the continuum formulation. In Section 4.5, we discuss shape DSA, in which design parameterization and design velocity field computations are first discussed. Shape DSA methods are followed. In Section 4.6, we discuss topology optimization, which is effective for generating structural layouts in support of concept design. We offer two case studies in Section 4.7 to showcase some of the typical design applications in practice. We assume that the performance measures cannot be expressed explicitly in design. Therefore, although simple problems are employed for illustration, readers should be mindful that in practice function evaluations are carried out numerically, such as using FEA.
Gradient-based algorithms have a solid mathematical background, in that Karush–Kuhn–Tucker (KKT) conditions are necessary for local minimal solutions. Under certain circumstances (for example, if the objective function is convex defined on a convex set), they can also be sufficient conditions. However, solving KKT conditions directly is usually very cumbersome, in that the equations are nonlinear, so practical algorithms aim to decrease the objective function value step by step instead (Leng, 2015). For unconstrained optimization, line search and trust region methods are most widely used. For more complicated constrained optimization, SQP, penalty, and projection methods can be utilized (Leng, 2015). Duality theory can also be a powerful tool by reformulating the primal problem in the dual space; the correspondent dual function is always concave, and sometimes the dual problem is much easier to solve. Nevertheless, all these methods require the gradient of a certain function with respect to its variables (Leng, 2015). An analytical form of the gradient is not guaranteed to exist, or could be very complicated. On the other hand, finite difference approximation of the gradient can be computationally costly for simulation-based optimization due to the increased number of objective function evaluations, which is usually characteristic of structural optimizations (Leng, 2015). For theoretical background in functional analysis, see Luenberger (1969). For numerical optimization theory and description of algorithms, see Nocedal and Wright (2006). Luenberger and Ye (2008) is also a recognized reference on mathematical programming theory.
SD is a simple, gradient-based optimizer that uses only first-order derivatives at the current design point to guide search. The well-established theoretical background of SD is quite intuitive: negative gradient is the SD direction in the neighborhood of the current point. SD is widely used on its own, and can also be integrated into other algorithms. The iterative scheme of SD with design variables updated at an iteration k is as follows:
[6.2]
where xk is the vector of design variables at iteration k, dk is the vector of design variable change direction, and scalar αk is a step length control parameter used to ensure improvement in the objective function. For SD, the step direction is the negative of the gradient of the function to be optimized.
[6.3]
The gradient , formed by partial derivatives of f with respect to each component of x, can be evaluated analytically or approximated by finite difference. When the vector norm of is zero, the necessary condition of a local minimum is satisfied and the algorithm has converged.
Given a step direction, a line search algorithm is used to identify the step size αk that produces a maximum reduction of the objective function. This is equivalent to finding the optimal value of scalar α for a function of only one variable, defined as
[6.4]
Various forms of line search have been developed; see Nocedal and Wright (2006) and Luenberger and Ye (2008) for additional details on line search algorithms.
Trust region method is another iterative algorithm effective for unconstrained numerical optimization. The terms “trust region” relates to a quadratic model mk of the objective function near the design point xk and the radius of that neighborhood Δk, developed using the Taylor-series expansion:
[6.5]
where , , and t is some scalar in the interval (0, 1). In the quadratic model function mk, an approximation Bk of the Hessian matrix (second-order partial derivatives of the objective function) is adopted.
The model function mk is minimized locally within the neighborhood, characterized by its radius Δk, ie, solving the subproblem:
[6.6]
Once a step pk is available, the performance of the quadratic model mk is evaluated using the ratio ρk as the predicted reduction of the objective function over the actual reduction of the model function:
[6.7]
If ρk is close to 1, indicating a good agreement between the model mk and the objective f, it is safe to expand the trust region for the next iteration. If ρk is positive but significantly smaller than 1, the radius of the trust region will remain the same; if ρk is close to zero or negative, the trust region will shrink. Also, the design point will be updated as if ρk is large enough.
The subproblem (Eq. [6.6]) is a constrained optimization problem, but it has some unique features that make it easier to solve. A straightforward method that is related to the idea of SD is to find pk along the direction of , called the Cauchy point. More detailed discussion on the trust region method is beyond the scope of this chapter but can be found in Nocedal and Wright (2006).
The trust region method has shown the power of using the quadratic model function in the optimization of general nonlinear objective functions. For general nonlinear programming problems with nonlinear objective functions and nonlinear equality and inequality constraints, the logic of the SQP method is inductive: at iteration k, an approximated quadratic programming is proposed as:
[6.8]
where is the matrix of second-order partial derivatives of the Lagrangian with respect to x; equality (index set ) or inequality (index set ) constraint is of the generic form. Note that in Eq. [6.8] the constraints are linearized. As in the trust region method, the can be approximated by a symmetric matrix Bk.
The solution of Eq. [6.8] can be obtained by the active-set or interior-point methods, which cannot be further explained due to space limitations. The new iterate at k + 1 is given by , where pk and λk+1 are the solution and Lagrange multiplier of Eq. [6.8]. The solution of Eq. [6.8] can further be modified using linear search technique as , since Eq. [6.8] is an approximated problem for pk. Methods of a trust region type have also been developed to solve Eq. [6.8].
Solving constrained mathematical programming problems is a challenging and fascinating task that is still under development. Other well-documented methods, like penalty and augmented Lagrangian methods and interior-point methods, can be found in Nocedal and Wright (2006) and Luenberger and Ye (2008).
The search for an optimal state is one of the most fundamental principles in our world. Optimization is to find the best solution to a certain designated problem. Numerical methods for optimizing the performance of engineering problems have been studied for many years. For optimization of an objective, different categories of optimization techniques namely “Gradient based” methods and “Non-gradient based” methods can be used. The basic idea behind the gradient based methods is that a given function reaches its extremes (minimum or maximum) in the direction of its gradient. Gradient based methods are in general computationally faster (they require fewer objective function evaluations in case of problems with low number of design variables) than non-gradient based methods. Their main drawback is that they might converge to local minima, and their convergence to global minima depends mainly on the chosen starting point (initial guess) by the user (Fig. 7a). On the other hand, non-gradient based methods, do not rely, most of the time, on strong mathematical basis and make use of neither the gradient or the second derivative of the objective function as a direction of descent. These methods work based on function evaluations alone. In principle, they attempt to mimic nature in order to find the optimum of the objective function. One of the key features of these algorithms is that they search from multiple points in the design space, instead of moving from a single point like gradient based methods. Thus, although in general there is no proof that these methods converge to global optima, experience entails that they converge to global optima in most cases (Fig. 7b).
Fig. 7. Optimization methods: (a) Gradient based (b) Non-gradient based.
The main drawback of non-gradient based methods is that these algorithms are generally computationally slower than the gradient based ones. Indeed, the non-gradient based algorithms can require thousands of function evaluations and, in some cases, become non-practical. In order to overcome these difficulties, the so-called Hybrid algorithms, which take advantage of the robustness of the non-gradient based methods and the fast convergence of the gradient based methods, have been proposed by different scholars. A set of analytically formulated rules and switching criteria can be coded into the program to automatically switch back and forth among the different algorithms as the iterative process advances. Each technique provides a unique approach with varying degrees of convergence, reliability and robustness at different stages during the iterative optimization process.
It is also interesting to introduce to the reader the multi-objective optimization problems vs. single-objective optimization problems. In a large number of problems, there exists a need to find optimal solutions due to more than one objective. In this case, a multi-objective approach must be employed. Engineering problems demanding low cost, high performance and low losses are an example of applications where this approach is needed. For the single-objective optimization problems, a unique optimal solution exists. However, for multi-objective optimization problems, there exist a set of compromised solutions, known as the Pareto-optimal solutions or non-dominated solutions, which are based on the competing objectives. The goal of multi-objective optimization is to find such set of solutions. Once the solutions are obtained, the designer can choose a final design with further considerations. Non-gradient based methods are more suitable for multi-objective optimization problems since they are able to find the entire range of Pareto-optimal solutions. Gradient based methods typically use a weighting method to combine different objectives into one for handling multi-objective optimization problems. In this case, results would be dependent on the chosen weighting coefficients.
An overview of different gradient based and non-gradient based optimization algorithms are provided in Appendix A. Clearly, knowledge about the nature of the problem is a requirement for choosing the most suitable optimization tool for an application.
3.2 Aerodynamic shape optimization of tall buildings using CFD
The goal of aerodynamic shape optimization is to accurately and efficiently determine surface shapes that attain optimal aerodynamic performance [1]. While the effects of geometric modifications to the shape of tall buildings, such as utilizing recessed or chamfered corners, etc. (discussed in Section 2.1) can significantly improve the aerodynamic response of tall buildings, a systematic approach for taking full advantage of aerodynamic shape optimization for buildings is not fully explored yet.
Experimental method using wind tunnel testing, provide the basis of the traditional “cut and try” approach for the design of new aerodynamic shapes. In this approach, several configurations are investigated in a wind tunnel and the one that yields the best aerodynamic performance is identified. The reason behind using wind tunnel testing for designing the new aerodynamic shapes is that the relation between the external shape of a building and the resulting intensity of the wind excitations is complicated and the improvements in wind effects that can be obtained by specific geometric modifications is difficult to predict without performing experiments [63]. Merrick and Bitsuamlak [64] examined the effects of building shape on wind loading patterns for high-rise buildings by analyzing various buildings with foot prints of square, circular, triangular, rectangular and elliptical shapes using wind tunnel database. Results of this research outlined the general wind loading characteristics of simple building shapes. The sensitivity of each building shape to vortex-shedding was determined. While it is possible to use wind tunnel testing for aerodynamic shape optimization, this approach is highly demanding due to the time and cost limitations for performing each test. As a matter of fact, only a limited number of possible configurations which are chosen based on engineering experience and judgement can be examined.
The use of computational simulation to scan many alternative designs has proved extremely valuable in practice. With the advances in computational fluid dynamics and computing power of modern computers, CFD has contributed to cut aerodynamic design cost and time scales by reducing the number of required wind tunnel tests. Currently, CFD is mainly used for estimating aerodynamic performance of a given structure configuration and it still suffers the limitation that it does not guarantee the identification of the best possible design. To ensure the realization of the true best design, the ultimate goal of computational simulation methods should not just be the analysis of prescribed shapes, but the automatic determination of the true optimum shape for the intended application. This is the underlying motivation for the combination of computational fluid dynamics with numerical optimization methods for aerodynamic shape optimization problems. Some of the earliest studies of such an approach were made by Hicks et al [65] and Hicks and Henne [66] for aircraft wing design. Aerodynamic shape optimization using CFD has been used for many years in aerospace [67,68] and automotive [69,70] industries and has recently become the subject of increasing interest in civil structures especially for aerodynamic design of the shape of tall buildings [71].
A general approach for aerodynamic shape optimization using CFD is shown in Fig. 8. The essential components of an aerodynamic shape optimization problem are objective functions, constraints, design variables that define the possible geometries, a flow solver, and a numerical optimization method. Examples of objective functions for aerodynamic shape optimization of tall buildings can be reducing the drag force and/or vortex-induced vibrations. The optimization algorithm finds the values of the geometric parameters that optimize the objective function while satisfying the constraints. These components should be selected carefully as they have a direct impact on the accuracy and efficiency of solution [1]. Samareh [72] provides summaries of shape parameterization techniques that can be used to define design variables. In general, it is important that the selected parameterization technique provides sufficient flexibility in order to realize truly optimal designs. It is also desirable that the number of parameters necessary to define the shape to be small so that a reasonable convergence rate of the optimization can be obtained.
As discussed in Section 3.1, various numerical optimization methodologies can be used for aerodynamic shape optimization applications. Hicks et al [65] were the first to apply gradient based methods to aerodynamic shape optimization problems. They used the method of feasible directions, which is based on conjugate gradients, to optimize airfoil shapes in transonic flow governed by the small-disturbance equation. Since this pioneering work, the application of gradient based methods to aerodynamic shape optimization problems remained an active area of research. However, as the aerodynamic shape optimization problem is a complex one with possibly many local optima, non-gradient based methods (also called Evolutionary Algorithms (EA)) such as Genetic Algorithm (GA) are more suitable for this application to ensure reaching the global optimum. Evolutionary Algorithms have the advantages such as robustness, suitability to parallel computing and simplicity in coupling CFD codes. Owing to these advantages over the non-gradient based methods, EAs have become increasingly popular in a broad class of design problems [73]. However, the implementation of aerodynamic shape optimization using CFD is intrinsically difficult for bluff civil structures due mainly to the turbulent flow field in which the structures are immersed, the high Reynolds number values, and the multi-objective nature of the design problem [74]. So, while there are clear advantageous for using non-gradient based optimization theories for aerodynamic shape optimization applications, the use of these algorithms for shape optimization of civil structures is often impractical since the objective functions and constraints are evaluated using CFD and an extremely high number of function calls which can easily reach the order of thousands, is required for obtaining the optimum solution.
One way of tackling the efficiency issue of evolutionary search methods is to use a Surrogate-Based Optimization (SBO) methodology in CPU-intensive aerodynamic shape optimization applications. The main idea in SBO is to utilize data sampling and surface fitting strategies to parameterize the space of possible solutions via a simple, computationally inexpensive model to be used for the purpose of numerical optimization. So, the whole optimization process is managed by the surrogate model outputs. This is often referred to as optimizing the response surface of the system. The basic process for a surrogate-based optimization consists of the following steps [75]:
1.
Defining a sampling plan for the design space: This contains both the samples required for constructing the surrogate model and some additional samples needed for verifying the surrogate model.
2.
Performing numerical simulations at design points selected from the sampling plan.
3.
Constructing the surrogate model: Different surrogate models can be used such as polynomial regression, Kriging, radial basis functions, neural networks, and support vector regression [76].
4.
Validating the model: This is to find out the predictive capabilities of the surrogate model.
5.
Updating the model: Repeating the last four steps until the desired model validation accuracy is obtained.
6.
Optimization: The search for the optimum is carried out on the surrogate model.
The approximation efficiency and its accuracy are major issues in SBO. If the problem has a high number of design variables, the construction of surrogate model may cause extremely high computational cost, which makes the approximation inefficient. Design of Experiment (DoE) can be used to reduce the number of design points [77]. However, by using this approach the global optimum might be missed due to the uncertainty at the predicted point which may mislead the optimization process in a wrong way. As a matter of fact, the selected DoE strategy is of paramount importance for achieving a satisfactory accuracy of the surrogate model. Generally, sampling plans that more evenly fill the design space can reduce the bias error of the surrogate. A plenty of surrogate methods, search algorithms and updating algorithms have been proposed in the literature for a variety of applications including CFD-based aerodynamic shape optimization [67,74,76,78].
It should also be noted that the repeated evaluation of the objective functions required during the optimization demands fast flow solutions [1]. As a matter of fact, the flow solver has a significant influence on the efficiency of the optimization. Buildings are immersed in very turbulent flows due to ground roughness effect which after interacting with buildings represent very wide temporal and spatial scales governed by highly non-linear Navier-Stokes differential equations which are computationally expensive to solve [79]. Detailed reviews of effective flow solvers for the Navier-Stokes equations are given by Pueyo [80]. Reynolds-Averaged Navier-Stokes (RANS) models can be used for obtaining fast and reasonably reliable simulations even for complex shapes and high Reynolds numbers, and can be used for obtaining both steady and unsteady aerodynamic measures (URANS). On the other hand, Large Eddie Simulation (LES) methods and/or hybrid RANS/LES methods are more accurate than RANS models for the prediction of the unsteady forces at the cost of higher computational burden. Note that the Navier-Stokes equations should be solved repeatedly for each updated geometry as the optimization progresses. This issue makes the problem further complicated. Therefore, an efficient approach can be to use URANS models for narrowing down the solution domain and then use the higher accuracy LES models when a greater level of detail is needed for example when comparing two final competing shapes [63].
One challenge in the application of numerical optimization to a building aerodynamic shape design problem is regeneration of the mesh for geometry variations from the initial base geometry. So, one very important feature of the CFD-based shape optimization technique should be to have the possibility of updating the CFD meshes automatically. Without this ability, the CFD model may need to be reprocessed before a suitable mesh with the required quality can be obtained. Also, the boundary conditions may need to be reapplied before the analysis can be launched. Mesh morphing techniques allow performing mesh modifications without changing the mesh topology and avoid a costly recalculation of the physics equations with every new CFD model. In these methods, the mesh is first imported into to the mesh morpher tool. The mesh is parameterized and then based on design constraints and realistic assumptions, deformation vectors are added to the model. After each optimization iteration the mesh morpher modifies the mesh with respect to the original topology, thus allowing a more robust and faster design optimization environment using CFD. Quality controls for meshing can be set to guarantee the quality of the CFD simulation results. This ensures that the algorithm is independent from external inputs during the optimization and makes mesh generation process faster. Mesh morphers have been used in many applications of shape optimization using CFD [81–83] mainly for aerospace and automotive structures. Wei et al [84] and Wei et al [85] developed a mesh morphing algorithm based on a Laplacian smoothing approach which was very robust even for marked geometric modifications which are expected in aerodynamic shape optimization of tall buildings.
Civil engineering structures which often have square or rectangular shapes provide a wide variety of applications for bluff body aerodynamics. The optimization of a trapezoidal bluff body was conducted by Burman et al [86]. This study focused on the response surface optimization of Navier-Stokes flow over a bluff body. The optimization objectives were to minimize the mean drag coefficient and to maximize the measure of mixing defined as the time average maximum negative velocity. Response surface methodology was utilized in the optimization process. Mack et al [87] used multiple surrogates for the shape optimization of a 2D trapezoidal bluff body model. Polynomial response surface and radial basis neural networks were used as surrogates. Objective functions were minimizing the total pressure loss coefficient and maximizing the mixing capability of time dependent flows over a 2-D bluff body.
Very recently, attempts have been made to use aerodynamic shape optimization in the field of tall building design. Bobby et al [71] proposed a framework for the aerodynamic shape optimization of tall buildings using CFD based on the concept of low-dimensional models for describing the global aerodynamic performance of tall buildings. A low-dimensional model was introduced which allows an efficient way for reducing the computational demand of the CFD simulations, thus expediting the optimization. A mesh morphing tool was used for automatically updating the mesh. The objective function was to minimize the drag force. Two-dimensional CFD simulation was used for each optimization iteration using Reynolds-Averaged Navier-Stokes (RANS) model. The low-dimensional model allowed extension of the results to the overall aerodynamic performance of tall buildings. Elshaer et al [88] proposed an optimization framework for aerodynamic shape optimization of tall buildings. The method couples the optimization algorithm using genetic algorithm (GA), computational fluid dynamics (CFD) solver, and the neural networks (NN) model. The objective was to reduce the drag force acting on a tall building by changing the shape of its corners. Large eddy simulation (LES) models were used for numerical simulation of the wind behavior. Bernardini et al [74] performed multi-objective aerodynamic shape optimization of tall buildings using a surrogate-based optimization method. The multiple objectives were minimizing the mean drag coefficient as well as the standard deviation of the lift coefficient. Ordinary Kriging surrogate model was used. A specifically developed strategy was adopted to update the Kriging models in order to perform additional CFD runs efficiently. Shell scripting, parallelized computations and mesh morphing algorithms were used to improve the efficiency and consistency of the framework.
These researches showed that in spite of recent advancements, there are still many challenges for the CFD-based aerodynamic shape optimization of tall buildings. These limited studies mainly focused on optimizing the dimensions of some predefined local corner modifications in tall buildings. To the author’s knowledge, there is no reported work till date on global shape optimization of tall buildings against wind load effects. Although this field of research is fairly new and is constantly improving, it is believed that full utilization of aerodynamic shape optimization could enable design of tall buildings that withstand wind load effects more efficiently. With the rapidly increasing use of the distributed parallel computing and adequate computing hardware, aerodynamic shape optimization is becoming more affordable than before.
The optimization algorithms can be divided into gradient-based algorithms and heuristic algorithms. The gradient is a vector composed of the partial derivatives of the objective function for each variable. The function grows fastest along the direction of the gradient and decreases fastest along the opposite direction of the gradient. Since the objective function and constraint function are not explicit functions of variables generally, it is difficult to solve the gradient. Therefore, some approximate methods are used. Lan et al. have done a lot of works, in which the gradient is approximated by finite difference method and realized by function of MATLAB Optimization Toolbox [28,31,33,38,86], which is a local optimization tool to find the minimum value. Liu et al. employed the adjoint method to approximate the gradient for the topology optimization by the method of Moving Asymptotes, which makes a series of subproblems with convexity and separable variables to approximate the nonlinear optimization problem [81,82]. Hence, the method of Moving Asymptotes may also be a feasible multi-parameter optimization algorithm.
The heuristic algorithms generally adopts stochastic iterative strategies, such as “genetic”, “crossover” and “mutation”, to screen according to objective function and constraint function, which narrows the solution space and makes an accurate local search. The Genetic Algorithm has been widely applied in multi-parameter optimization [48,49,92]. Most of them are implemented in MATLAB Global Optimization Toolbox. Moreover, there are some achievements in the optimization and variation of Genetic Algorithm. Pham, Wang, et al. employed the Genetic Algorithm with non-dominated sorting program to derive the load–deflection relation for better efficiency, which is realized by simplifying multiple objectives into a virtual fitness function [37,53,76,94]. The Particle Swarm Optimization Algorithm captures the optimal solution through the cooperation and information sharing among individuals in particle swarm without many parameters to adjust, which can be adopted to optimize the structural parameters of CCFM [99]. Each particle searches the optimal solution individually in the design domain, that is, individual extreme value, and shares it with other particles in the particle swarm. The optimal value of individual extreme value is regarded as the global optimal solution.
Different from the uniform section beam in all the above CCFMs, the interpolation optimization generate the flexible beams with variable section to generate constant force. So far, only cubic spline interpolation (see [42,54,93]) is an available case. The neutral axis of cubic spline curve based on the five interpolation points , , , and is shown in Fig. 18. The radius of each section perpendicular to the neutral axis at each interpolation point is not all the same value. The variation of section radius from one interpolation point to adjacent interpolation points is continuous, linear and monotonic. The synthesis of CCFM by interpolation optimization is realized by MATLAB Global Optimization Toolbox. Because this type curved beam has some rapidly changing angles, such as the angle at , much attention should be paid to avoid the stress concentration caused by improper section radius.
Fig. 18. Cubic spline interpolation with five different radius of the interpolation circles, reproduced from [93].
The calculation of derivatives or the Hessian matrix is the key to solving optimization problems using gradient-based algorithms. Researchers have used various calculation methods including numerical perturbation, sensitivity equation, and adjoint method, as summarized in Table 7.11.
Table 7.11. Comparison of calculation methods for gradients and Hessian matrix.
Methods
Calculation principles
Characteristics
Numerical perturbation
Small perturbations of the model parameters and calculation of the production responses
Easy to implement; expensive computational cost; unsuitable for large-scale optimization problems
Sensitivity equation
Differentiation of the flow and transport equations
Difficult to obtain analytical expressions for nonlinear optimization problems
Adjoint method
Optimal control theories and calculus of variations
Easy to implement; dependent on reservoir simulators; hard to transplant elsewhere
From Hou, J., Zhou, K., Zhang, X.S., Kang, X.D., Xie, H., 2015. A review of closed-loop reservoir management. Petroleum Science 12,114–128. https://ds.doi.org/10.1007/s12182-014-0005-6.
Generally, gradient-based algorithms can be classified into two categories. One is the first-order methods, which only require the derivative information. The other category is the second-order methods, which not only require the derivative information but also require the Hessian matrix. Representative methods include Gauss–Newton, Levenberg–Marquardt, sequential quadratic programming (Barnes et al., 2007), and the limited memory Broyden Fletcher Goldfarb Shanno method (LBFGS).
When using the Gauss–Newton method for automatic history matching problems, Wu et al. (1999) introduced an artificially high variance of measurement errors at early iterations to damp the changes in model parameters and thus avoid undershooting or overshooting. Tan and Kalogerakis (1991) pointed out that the standard Gauss–Newton and Levenberg–Marquardt methods require the computation of all sensitivity coefficients in order to formulate the Hessian matrix, which seems impossible in reality due to the large number of unknown parameters relative to limited available measurements. In order to eliminate this problem, the quasi-Newton method was introduced by researchers. This method only requires the gradient of the objective function which can be computed from a single adjointsolution as done in Zhang et al. (2002). In order to further improve the computational efficiency and robustness of the LBFGS method, Gao and Reynolds (2006) proposed a new line search strategy, rescaled the model parameters, and applied damping factors to the production data. They also noticed that the new line search strategy had to satisfy the strong Wolfe conditions at each iteration, or the convergence rate would decrease significantly.
The Karhunen–Loeve expansion can create a differentiable parameterization of the numerical model in terms of a small set of independent random variables and deterministic eigenfunctions. With this expansion, the gradient-based algorithms can be applied while honoring the two-point statistics of the geological models (Gavalas et al., 1976). In order to further extend the existing gradient-based history matching techniques to deal with complex geological models characterized by multiple-point geostatistics, Sarma et al. (2007) applied a kernel principal component analysis method to model permeability fields. This method can preserve arbitrary high-order statistics of random fields, and it is able to reproduce complex geology while retaining reasonable computational requirements.