Methods for Parametric Design Optimization

ANSYS, Inc. Release 2025 R2

Southpointe July 2025
2600 Ansys Drive
Canonsburg, PA 15317 ANSYS, Inc. and
ansysinfo@[Link] ANSYS Europe,
Ltd. are UL
[Link] registered ISO
(T) 724-746-3304 9001: 2015
(F) 724-514-9494 companies.
Copyright and Trademark Information

© 2025 ANSYS, Inc. Unauthorized use, distribution or duplication is prohibited.

ANSYS, Ansys Workbench, AUTODYN, CFX, FLUENT and any and all ANSYS, Inc. brand, product, service and feature
names, logos and slogans are registered trademarks or trademarks of ANSYS, Inc. or its subsidiaries located in the
United States or other countries. ICEM CFD is a trademark used by ANSYS, Inc. under license. CFX is a trademark
of Sony Corporation in Japan. All other brand, product, service and feature names or trademarks are the property
of their respective owners. FLEXlm and FLEXnet are trademarks of Flexera Software LLC.

Disclaimer Notice

THIS ANSYS SOFTWARE PRODUCT AND PROGRAM DOCUMENTATION INCLUDE TRADE SECRETS AND ARE CONFID-
ENTIAL AND PROPRIETARY PRODUCTS OF ANSYS, INC., ITS SUBSIDIARIES, OR LICENSORS. The software products
and documentation are furnished by ANSYS, Inc., its subsidiaries, or affiliates under a software license agreement
that contains provisions concerning non-disclosure, copying, length and nature of use, compliance with exporting
laws, warranties, disclaimers, limitations of liability, and remedies, and other provisions. The software products
and documentation may be used, disclosed, transferred, or copied only in accordance with the terms and conditions
of that software license agreement.

ANSYS, Inc. and ANSYS Europe, Ltd. are UL registered ISO 9001: 2015 companies.

U.S. Government Rights

For U.S. Government users, except as specifically granted by the ANSYS, Inc. software license agreement, the use,
duplication, or disclosure by the United States Government is subject to restrictions stated in the ANSYS, Inc.
software license agreement and FAR 12.212 (for non-DOD licenses).

Third-Party Software

See the legal information in the product help files for the complete Legal Notice for ANSYS proprietary software
and third-party software. If you are unable to access the Legal Notice, contact ANSYS, Inc.

Published in the U.S.A.

Table of Contents
1. Introduction ............................................................................................................................................ 1
2. Design of Experiments ............................................................................................................................ 3
2.1. Deterministic DoE schemes ............................................................................................................... 3
2.1.1. Full factorial design .................................................................................................................. 4
2.1.2. Star points ............................................................................................................................... 4
2.1.3. Central composite design ......................................................................................................... 5
2.1.4. Box-Behnken design ................................................................................................................ 6
2.1.5. Koshal designs ......................................................................................................................... 7
2.1.6. D-optimal designs .................................................................................................................... 8
2.1.7. Orthogonal Arrays .................................................................................................................... 9
2.2. Random and quasi-random DoE schemes ....................................................................................... 11
2.2.1. Monte Carlo Simulation .......................................................................................................... 11
2.2.2. Latin Hypercube Sampling ..................................................................................................... 12
[Link]. Correlation optimized Latin Hypercube Sampling .......................................................... 14
[Link]. Space-filling Latin Hypercube Sampling ......................................................................... 14
2.2.3. Space-filling designs .............................................................................................................. 15
2.2.4. Sobol Sequences .................................................................................................................... 17
2.3. Deterministic vs. random DoE schemes ........................................................................................... 18
2.4. Glossary ......................................................................................................................................... 19
3. Metamodeling techniques .................................................................................................................... 21
3.1. Regression models for scalar outputs .............................................................................................. 21
3.1.1. Polynomial regression ............................................................................................................ 22
3.1.2. Box-Cox transformation .......................................................................................................... 23
3.1.3. Moving Least Squares approximation ..................................................................................... 24
3.1.4. Kriging ................................................................................................................................... 26
3.1.5. Radial basis functions ............................................................................................................. 27
3.1.6. Support Vector Regression ..................................................................................................... 30
3.1.7. Sparse Grid Approximation ..................................................................................................... 32
3.1.8. Genetic Aggregation Response Surface .................................................................................. 37
3.1.9. Neural networks ..................................................................................................................... 40
[Link]. Feedforward neural networks ........................................................................................ 43
[Link]. Radial basis function networks ....................................................................................... 46
[Link]. Deep Feed Forward Network ......................................................................................... 50
[Link].1. Hyperparameter search ......................................................................................... 51
[Link].2. Model training ...................................................................................................... 51
[Link].3. Network layout ..................................................................................................... 52
3.1.10. Deep Infinite Mixture Gaussian Process (DIM-GP) .................................................................. 53
3.2. Model quality measures .................................................................................................................. 54
3.2.1. Residual analysis .................................................................................................................... 55
[Link]. Residual sum of squares ................................................................................................. 55
[Link]. Mean squared error (MSE) .............................................................................................. 55
[Link]. Root mean squared error (RMSE) .................................................................................... 55
[Link]. Relative root mean squared error ................................................................................... 55
[Link]. Maximum residual ......................................................................................................... 56
[Link]. Maximum Relative Residual ........................................................................................... 56
[Link]. Relative Maximum Absolute Error .................................................................................. 56
[Link]. Relative Average Absolute Error ..................................................................................... 56
[Link]. PRESS residuals for polynomial regression ...................................................................... 57
[Link]. R2 for prediction of polynomial regression .................................................................... 57

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. iii
Methods for Parametric Design Optimization

3.2.2. Generalized cross-validation and Akaike"s final prediction error .............................................. 58

3.2.3. k-fold cross-validation ............................................................................................................ 58
3.2.4. Coefficient of Determination (R2) ............................................................................................ 59
3.2.5. Coefficient of Prognosis .......................................................................................................... 61
3.3. Model Sensitivity measures ............................................................................................................. 63
3.3.1. Variance based sensitivity analysis .......................................................................................... 63
[Link]. ANOVA .......................................................................................................................... 63
[Link].1. The confidence interval of the regression coefficients ............................................ 64
[Link].2. The significance of a regression coefficient ........................................................ 65
[Link]. Principal Component Analysis ........................................................................................ 66
[Link]. Coefficient of Correlation ............................................................................................... 67
[Link]. Coefficient of Importance .............................................................................................. 67
3.4. MOP principle and competition ...................................................................................................... 68
3.4.1. MOP residual plot ................................................................................................................... 70
3.4.2. Local MOP error measures ...................................................................................................... 71
3.5. Adaptive metamodeling ................................................................................................................. 73
3.5.1. Adaptive Metamodel of Optimal Prognosis (AMOP) ................................................................ 73
3.5.2. Genetic Aggregation Response Surface Refinement ................................................................ 75
3.6. Metamodels Glossary ...................................................................................................................... 77
4. Optimization methods .......................................................................................................................... 81
4.1. Optimization setup ......................................................................................................................... 81
4.2. Single-objective optimization ......................................................................................................... 84
4.2.1. Accompanying example: optimization of a damped oscillator ................................................. 84
4.2.2. Gradient-based methods ........................................................................................................ 86
[Link]. Non-Linear Programming by Quadratic Lagrangian (NLPQLP) ......................................... 87
[Link]. Mixed-Integer Sequantial Quadratic Programming (MISQP) ............................................ 89
[Link]. Leapfrog optimizer for constrained minimization (LFOPC) .............................................. 90
[Link]. DAKOTA Coliny Solis-Wets .............................................................................................. 91
[Link]. DAKOTA CONMIN .......................................................................................................... 91
[Link]. DAKOTA OPT++ Finite Differences Newton (FDN) ........................................................... 91
[Link]. DAKOTA OPT++ Quasi Newton (QN) ............................................................................... 91
[Link]. DAKOTA OPT++ Polak-Ribiere Conjugate Gradient (PR) ................................................... 92
4.2.3. Pattern search methods .......................................................................................................... 92
[Link]. Downhill simplex ........................................................................................................... 92
[Link]. Hooke-Jeeves Pattern Search ......................................................................................... 93
[Link]. DAKOTA OPT++ Parallel Direct Search (PDS) ................................................................... 94
[Link]. DAKOTA Asynchronous Parallel Pattern Search (APPS) .................................................... 94
[Link]. DAKOTA Coliny DIRECT .................................................................................................. 94
[Link]. DAKOTA NCSU DIRECT ................................................................................................... 95
[Link]. DAKOTA Coliny Pattern Search ....................................................................................... 96
4.2.4. Response surface based methods ........................................................................................... 96
[Link]. Adaptive Response Surface Method (ARSM) ................................................................... 97
[Link]. Optimization using the Metamodel of Optimal Prognosis ............................................... 99
[Link]. Adaptive Metamodel of Optimal Prognosis (AMOP) ...................................................... 100
[Link]. Adaptive Single-Objective ........................................................................................... 101
[Link]. Efficient Global Optimization ....................................................................................... 104
[Link]. Probabilistic Inference for Bayesian Optimization ......................................................... 107
4.2.5. Nature-inspired methods ...................................................................................................... 109
[Link]. Evolutionary algorithms ............................................................................................... 110
[Link]. Darwin algorithm ........................................................................................................ 113
[Link]. EVOLVE ........................................................................................................................ 114

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
iv of ANSYS, Inc. and its subsidiaries and affiliates.
 Methods for Parametric Design Optimization

[Link]. DAKOTA Coliny Evolutionary Algorithm (EA) ................................................................. 114

[Link]. Covariance Matrix Adaptation ...................................................................................... 115
[Link]. Particle Swarm Optimization ........................................................................................ 115
[Link]. Stochastic Design Improvement ................................................................................... 117
[Link]. Adaptive Simulated Annealing ..................................................................................... 118
[Link].1. Algorithm ........................................................................................................... 119
[Link].2. Acceptance function ........................................................................................... 119
[Link].3. Sampling algorithm ............................................................................................ 120
[Link].4. Cooling schedule ................................................................................................ 120
[Link].5. Stopping criterion ............................................................................................... 121
[Link].6. Re-annealing ...................................................................................................... 121
[Link].7. Some comments ................................................................................................. 122
[Link]. Differential Evolution ................................................................................................... 124
4.2.6. Hybrid methods ................................................................................................................... 127
[Link]. One-Click Optimization (OCO) ...................................................................................... 127
4.3. Multi-objective optimization ......................................................................................................... 128
4.3.1. Pareto optimization .............................................................................................................. 130
[Link]. Pareto optimality ......................................................................................................... 130
[Link]. Analysis of conflicting objectives .................................................................................. 132
4.3.2. Weighted sum ...................................................................................................................... 132
4.3.3. ε-Constraint ......................................................................................................................... 133
4.3.4. Response Surface based methods ......................................................................................... 133
[Link]. Adaptive Metamodel of Optimal Prognosis (AMOP) ...................................................... 133
[Link]. Adaptive Multiple-Objective ........................................................................................ 134
4.3.5. Natured inspired methods .................................................................................................... 136
[Link]. Evolutionary algorithms ............................................................................................... 136
[Link]. Darwin algorithm ........................................................................................................ 137
[Link]. Particle Swarm Optimization ........................................................................................ 138
[Link]. Non-dominated Sorting Genetic Algorithm II ............................................................... 138
[Link]. Multi-Objective Genetic Algorithms ............................................................................. 141
4.3.6. Hybrid methods ................................................................................................................... 142
[Link]. One-Click Optimization (OCO) ...................................................................................... 142
4.3.7. Performance metrics ............................................................................................................ 142
[Link]. Number of nondominated points ................................................................................. 142
[Link]. Spread ........................................................................................................................ 142
[Link]. Standard deviation of crowding distance ..................................................................... 142
[Link]. Min/Max of objectives ................................................................................................. 143
[Link]. Hypervolume .............................................................................................................. 143
[Link]. Number of common points .......................................................................................... 143
[Link]. Number of new nondominated solutions ..................................................................... 144
[Link]. Number of old dominated solutions n(Q) ..................................................................... 144
[Link]. Consolidation ratio ...................................................................................................... 144
[Link]. Improvement ratio ..................................................................................................... 144
4.4. Optimization Glossary ................................................................................................................... 145

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. v
 Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
vi of ANSYS, Inc. and its subsidiaries and affiliates.
List of Figures
2.1. Full factorial (left) and full combinatorial design (right) in a 3-dimensional space. ..................................... 4
2.2. Star points for in a 3-dimensional space. .......................................................................................... 5
2.3. Central composite design in a 3-dimensional space. ................................................................................ 6
2.4. Box and Behnken design in a 3-dimensional space. ................................................................................. 7
2.5. Linear Koshal design (left) and quadratic Koshal (right) in a 3-dimensional space. ..................................... 8
2.6. Linear D-optimal design (left) and quadratic D-optimal (right) in a 3-dimensional space. .......................... 8
2.7. Monte Carlo Simulation in a 2-dimensional space. ................................................................................. 11
2.8. Latin Hypercube Sampling in a 2-dimensional space. ............................................................................. 13
2.9. Advanced correlation optimized LHS (left) and space filling LHS using the centered L2-discrepancy ap-
proach (right) in a 2-dimensional space. ...................................................................................................... 14
2.10. Six space-filling designs: 5 points in a 2-dimensional space. .................................................................. 16
2.11. Sobol sequences in a 2-dimensional space with 63 samples (left)and eight dimensions with 127 samples
(right) ......................................................................................................................................................... 18
2.12. Dimension reduction for a nonlinear function with five inputs based on Latin Hypercube Sampling (left)
with 100 samples and full factorial design (right) with samples ....................................................... 19
3.1. Box-Cox transformation for different values of .................................................................................... 23
3.2. Local weighting of support point values of the MLS approximation ........................................................ 24
3.3. MLS approximation using a smoothing Gaussian kernel as weighting function (left) and an interpolating
kernel (right) depending on the influence radius ......................................................................................... 25
3.4. Smoothing Kriging approximation for noisy data ................................................................................... 27
3.5. Different Radial Basis Functions formulations depending on the radius .................................................. 28
3.6. Radial Basis Function interpolation for the different basis types. ............................................................. 29
3.7. Radial Basis Function extension for noisy data. ...................................................................................... 29
3.8. Regression line for a group of sample points with a tolerance of , which is characterized by slack variables
and . .................................................................................................................................................... 31
3.9. Piecewise linear hierarchical basis (from the level 0 to the level 3). .......................................................... 33
3.10. Interpolation with the Hierarchical Basis. ............................................................................................. 33
3.11. Tensor product approach to generate the piecewise bilinear basis functions ................................ 34
3.12. Tensor product of linear basis functions for a two-dimensional problem ............................................... 35
3.13. Discretization of piecewise linear basis , and ...................................................................... 36
3.14. Sparse Grid example of 2D response surface. ....................................................................................... 37
3.15. Cross-validation error at the second sample point: . ...................................................................... 39
3.16. Schematic of a neural network with 2 inputs and a hidden layer of 4 neurons with activation function
. ............................................................................................................................................................... 41
3.17. Neural Network examples ................................................................................................................... 42
3.18. Schematic view of a Deep Feed Forward Network with 2 inputs and 3 hidden layers of 4 neurons with
activation function f. ................................................................................................................................... 50
3.19. Schematic view on cross-validation procedure used for the training and evaluation of the DFFN mod-
el. ............................................................................................................................................................... 52
3.20. Schematic representation of DIM-GP ................................................................................................... 53
3.21. Subspace plot of the investigated nonlinear function (Equation 3.73) and convergence of the CoD
measures with increasing number of learning points ................................................................................... 60
3.22. Convergence of the CoP measure by using MLS approximation compared to the polynomial CoD
measure ..................................................................................................................................................... 62
3.23. Definition of ANOVA value .............................................................................................................. 64
3.24. ANOVA value with confidence interval ............................................................................ 65

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. vii
Methods for Parametric Design Optimization

3.25. CoP values of different input variable combinations and approximation methods obtained with the
analytical nonlinear function ....................................................................................................................... 68
3.26. - and - subspace plots of the MOP of the nonlinear analytical function given in Equation
3.73 ............................................................................................................................................................ 69
3.27. Residual plot of the MOP-post-processing: original vs. the approximated values (left) and the sample
CoP indicated for each design (right) ........................................................................................................... 71
3.28. MOP-post-processing: local CoP calculated as weighted averaging of the sample CoPs ......................... 72
3.29. Estimated support density of a design set with almost uniform distribution (left) and a pure random
distribution with regions with high and low density (right) .......................................................................... 72
3.30. Adaptive MOP: approximation history plot with the convergence of the global and minimum local CoP
values for each response ............................................................................................................................. 73
3.31. Adaptive MOP: Local refinement considering the local approximation quality ...................................... 74
3.32. Adaptive MOP: Local refinement considering constraints of the response values .................................. 74
3.33. Illustration of the UP distribution for a kriging surrogate (left). Dashed lines: CV sub-models predictions,
solid red line: master model prediction, horizontal bars: local UP distribution at = -1.8 and = 0.2, black
squares: design points. ................................................................................................................................ 75
3.34. Uncertainty quantification based on the UP distribution for a kriging surrogate. Blue solid line: master
model prediction , light blue area: region delimited by ................................................ 75
3.35. Absolute predicted error in a 1-dimensional problem. Circles: design points used to build the
metamodel, diamond: 1st refinement point ................................................................................................. 76
3.36. Absolute weighted predicted error in a 1-dimensional problem. Circles: design points used to build the
metamodel, diamond: 1st refinement point, triangle: 2nd refinement point .................................................. 76
4.1. Definition of optimization parameter properties in optiSLang ................................................................ 81
4.2. Definition of objective functions and optimization constraints in optiSLang ........................................... 82
4.3. Recommended flowchart for single-objective optimization ................................................................... 84
4.4. Damped oscillator: system properties and oscillation behavior ............................................................... 84
4.5. Objective function of the damped oscillator obtained by using an estimate of the maximum amplitude
from the envelope curve (left) and by using a coarse time discretization of the displacement curve (right) .... 86
4.6. Iterative search of a gradient-based method using a quadratic approximation of the objective func-
tion ............................................................................................................................................................ 87
4.7. Convergence of the NLPQL optimizer for the oscillator optimization problem by using the smooth ob-
jective function based on the envelope estimate (left) and by using the noisy objective function from the
time discretization (right) ............................................................................................................................ 88
4.8. Main steps of the downhill simplex algorithm in a 2-dimensional [Link] coefficients used in this ex-
ample are not fixed and can vary. ................................................................................................................ 93
4.9. Subdivision steps of DIRECT algorithm (DAKOTA 5.0 Reference Manual, p. 70) ......................................... 95
4.10. Global polynomial approximation of the maximum amplitude (left, ) and the damped eigen-
frequency (right, ) of the optimized oscillator using a quadratic basis and a full factorial design
scheme ...................................................................................................................................................... 97
4.11. Adaptation of the polynomial approximation scheme inside the Adaptive Response Surface Method .... 98
4.12. Convergence of the Adaptive Response Surface Method for the damped oscillator by using the noisy
objective function: designs used for the adaptation (left) and modification of the local DoE bounds during
the iteration (right) ..................................................................................................................................... 99
4.13. Approximation of the maximum amplitude (left, ) and of the damped eigen-frequency (right,
) by using the Metamodel of Optimal Prognosis with 100 Latin Hypercube samples ................... 100
4.14. ASO workflow ................................................................................................................................... 102
4.15. From the current OSF to the new OSF ................................................................................................ 103

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
viii of ANSYS, Inc. and its subsidiaries and affiliates.
 Methods for Parametric Design Optimization

4.16. Example of 1-dimensional EGO algorithm. Dashed red line: real function, dashed black line: value of
the current minimum, solid line: Kriging prediction ( ), light blue area: Kriging variance ( ), black squares:
design points. ........................................................................................................................................... 106
4.17. Beginning of Bayesian Optimization. ................................................................................................. 108
4.18. Illustration of the Bayesian optimization iterations and the associated acquisition function. ................ 108
4.19. General flowchart of nature-inspired optimization algorithms ............................................................ 110
4.20. Evolution strategies versus genetic algorithms ................................................................................... 111
4.21. Convergence of the evolutionary algorithm with global search (left) and local search (right) for the
damped oscillator with noisy objective function ........................................................................................ 113
4.22. Convergence of the covariance matrix adaptation for the rosenbrock function ................................... 115
4.23. Update of a particle position in the Particle Swarm Optimization using a combination of the old ve-
locity , the direction to the local best position and the direction to the global best position .......... 116
4.24. Update scheme of the Stochastic Design Improvement approach by generating a random sampling
around the best design of a previous iteration ........................................................................................... 118
4.25. OCO workflow ................................................................................................................................... 128
4.26. Design space (left) and objective space (right) ................................................................................... 129
4.27. Recommended flowchart for multi-objective optimization ................................................................. 130
4.28. Pareto dominating (filled circles) and dominated designs (unfilled circles) including the resulting Pareto
frontier of two conflicting objectives. ( and dominate but are indifferent to each other) ...................... 131
4.29. Conflicting objectives of the damped oscillator (maximum amplitude vs. eigen-frequency) including
cluster analysis of the Pareto optimal designs in the objective space (left) and in the design space (right) .... 132
4.30. Adaptive MOP: Pareto frontier update for the damped oscillator example using the default space-filling
update criterion (left) and by using the best compromise criterion (right) .................................................. 134
4.31. AMO workflow .................................................................................................................................. 135
4.32. Pareto frontier for the damped oscillator (left) and Pareto optimal designs plotted in the design space
(right) obtained by evolutionary algorithms .............................................................................................. 137
4.33. Elitist non-dominated sorting genetic algorithm (NSGA-II).The shaded blocks are not the part of original
NSGA-II but additions to avoid Pareto drift. ................................................................................................ 139
4.34. Elitism in NSGA-II. .............................................................................................................................. 140
4.35. Illustration of non-domination criterion, Pareto optimal set, and Pareto optimal front. ........................ 140
4.36. Hypervolume metric of two conflicting objectives ( and ) with as reference point. .................... 143

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. ix
 Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
x of ANSYS, Inc. and its subsidiaries and affiliates.
List of Tables
2.1. Required number of designs depending on the input dimension for a linear and quadratic approximation
using Koshal, D-optimal, full factorial, central composite (CCD) and Box-Behnken design schemes .................. 3
2.2. orthogonal array which consists of 8 designs in 7-dimensional unit cube ...................................... 9
2.3. Interaction table for the orthogonal array ................................................................................... 10
2.4. Description of space-filling algorithms .................................................................................................. 16
2.5. Recommended number of samples for the Sobol sequences depending on the input dimension ........... 18
3.1. Available regression models for scalar outputs and recommended application w.r.t. the number of
training data points .................................................................................................................................... 22
4.1. Available optimization algorithms with application for single and multi-objective optimization and
possible discrete inputs variables ................................................................................................................ 82

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. xi
 Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
xii of ANSYS, Inc. and its subsidiaries and affiliates.
Chapter 1: Introduction
Optimization and robustness analyses have become important tools for the virtual development of in-
dustrial products. In parametric optimization, the optimization variables are systematically modified by
mathematical algorithms in order to get an improvement of an existing design or to find a global op-
timum. The design variables (also named input parameters) are defined by their lower and upper bounds
or by several possible discrete values. In real world industrial optimization problems, the number of
design variables can often be very large. Unfortunately, the efficiency of mathematical optimization al-
gorithms decreases with increasing number of design variables. For this reason, several methods are
limited to a moderate number of variables, such as gradient based and Adaptive Response Surface
Methods. With the help of sensitivity analysis, the designer identifies the variables which contribute
most to a possible improvement of the optimization goal. Based on this identification, the number of
design variables may be dramatically reduced and an efficient optimization can be performed. Additional
to the information regarding important variables, sensitivity analysis may help to decide, if the optimiz-
ation problem is formulated appropriately and if the numerical CAE solver behaves as expected.

A modern approach to search for better designs or to compute the "best" design has to introduce all
available engineering know-how and has to automate a multidisciplinary optimization process. By
specifying the design criteria as objectives and constraints and specifying the space of all possible
designs with optimization parameters, a framework for numerical optimization can be defined. Part of
the challenge of defining a multidisciplinary optimization problem will be the communication of different
design groups about conflicting objectives, fixed criteria or weighted compromise objective functions.
Consequently, the degree of non-linearity has to be taken into account in the optimization process.
Because of that, the resulting optimization problem may become very noisy, very sensitive to design
changes or ill-conditioned for mathematical function analysis (non-differentiable, non-convex, non-
smooth).

That defines the requirements to a modern optimization tool. Arbitrary solvers have to be connected,
optimization strategies for smooth as well as for non-smooth or even ill-posed problems have to be
available.

Besides the problem to find an optimal design, the evaluation of the robustness, i.e. the sensitivity of
unavoidable scatter of design variables due to the structural response, becomes more and more import-
ant. Very often, optimized designs tend to be very sensitive to small (sometimes random) fluctuations
of parameters. Such phenomena may occur due to system instabilities like bifurcation problems in the
structure. Design robustness can be checked by applying a systematic perturbation analysis based on
a randomly generated design sample set. Statistics on the sample set allows the evaluation of robustness
of the optimized design.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 1
 Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
2 of ANSYS, Inc. and its subsidiaries and affiliates.
Chapter 2: Design of Experiments
The following topics are available:
2.1. Deterministic DoE schemes
2.2. Random and quasi-random DoE schemes
2.3. Deterministic vs. random DoE schemes
2.4. Glossary

2.1. Deterministic DoE schemes

In order to perform a global sensitivity analysis, the space of the design variables has to be scanned by
discrete realizations. Each realization is one set of values belonging to the specified inputs. As design
exploration often deterministic Designs of Experiments (DoE) are applied (Myers and Montgomery 2002).
These design schemes are mainly based on a regular arrangement of the designs, as in the full factorial
design, where the number of designs increases exponentially with increasing dimension. Fractional
factorial designs use only a part of the full factorial designs and are often optimized to fit linear or
quadratic approximation functions with minimum errors. In the following section, all available design
schemes are introduced. In Table 2.1 (p. 3) the required number of designs for the most common
deterministic schemes are shown depending on the input dimension for different DoE schemes.

Table 2.1: Required number of designs depending on the input dimension for a linear and quadratic
approximation using Koshal, D-optimal, full factorial, central composite (CCD) and Box-Behnken
design schemes

Linear approximation Quadratic approximation

Dimension Koshal D-optimal Factorial Koshal D-optimal Factorial CCD Box-Behnken
1 2 2 2 3 3 3 - -
2 3 4 4 6 9 9 9 -
3 4 6 8 10 15 27 15 13
4 5 8 16 15 23 81 25 25
5 6 9 32 21 32 243 43 41
6 7 11 64 28 42 729 77 49
7 8 12 128 36 54 2187 143 57
8 9 14 256 45 68 6561 273 113
9 10 15 512 55 83 19683 531 121
10 11 17 1024 66 99 59049 1045 161

Remark 1 (D-optimal design). In Table 2.1 (p. 3) the number of designs for the D-optimal design schemes
are given for the optiSLang implementation, in LS-OPT the implementation requires slightly less designs.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 3
Design of Experiments

References
Myers, R., and D. C. Montgomery. 2002. Response Surface Methodology. 2nd ed. John Wiley & Sons, Inc.

2.1.1. Full factorial design

Available in: LS-OPT, ModelCenter, optiSLang

Figure 2.1: Full factorial (left) and full combinatorial design (right) in a 3-dimensional space.

The full factorial design is a multi-dimensional grid, which consists of points (levels) in each di-
mension. The total number of samples is given as , where is the number of input parameters.

The full factorial design can represent linear and multi-linear terms. For more than 2 levels, also
quadratic terms can be considered.

The full combinatorial design is a special case for discrete parameters, where the continuous parameters
use the user specified number of levels and the discrete parameters consider all possible discrete
values.

2.1.2. Star points

Available in: optiSLang

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
4 of ANSYS, Inc. and its subsidiaries and affiliates.
 Deterministic DoE schemes

Figure 2.2: Star points for in a 3-dimensional space.

The star point design consists of a center point and points in each axis direction (in general
on each axis). For inputs parameters the number of designs reads .

2.1.3. Central composite design

Available in: DesignXplorer, LS-OPT, ModelCenter, optiSLang

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 5
Design of Experiments

Figure 2.3: Central composite design in a 3-dimensional space.

The central composite design is a joined design consisting of a full factorial design and a star
point design, which results in designs in total. The full factorial design is defined in
and bounds, whereas the star point design is scaled by a factor . As discussed in (Myers and
Montgomery 2002), different values are possible.

Remark 1. In optiSLang the inner star points of the central composite design is scaled with -factors
equal to one. In LS-OPT, is assumed.

References
Myers, R., and D. C. Montgomery. 2002. Response Surface Methodology. 2nd ed. John Wiley & Sons,
Inc.

2.1.4. Box-Behnken design

Available in: DesignXplorer, ModelCenter, optiSLang

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
6 of ANSYS, Inc. and its subsidiaries and affiliates.
 Deterministic DoE schemes

Figure 2.4: Box and Behnken design in a 3-dimensional space.

Brief description: The Box-Behnken design (Box and Behnken 1960 (p. 7)) is sufficient to fit a quad-
ratic model, which may contain squared terms, products of two factors, linear terms and an intercept.
The property of "missing corners" can be useful if these should be avoided due to physical or economic
constraints, because potential of data loss in those cases can be prevented.

Limitations: The design scheme requires at least three input parameters and results in
designs.

References
Box, George EP, and Donald W Behnken. 1960. "Some New Three Level Designs for the Study of
Quantitative Variables." Technometrics 2 (4): 455–75.

2.1.5. Koshal designs

Available in: LS-OPT, optiSLang

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 7
Design of Experiments

Figure 2.5: Linear Koshal design (left) and quadratic Koshal (right) in a 3-dimensional space.

Koshal designs are fractional factorial designs which are saturated to model exactly a linear or quad-
ratic response. The number of designs is in case of a linear function and
for a quadratic function.

2.1.6. D-optimal designs

Available in: LS-OPT, optiSLang

Figure 2.6: Linear D-optimal design (left) and quadratic D-optimal (right) in a 3-dimensional
space.

The D-optimal schemes are generated as a subset of designs of the (linear) and (quadratic) full
factorial schemes. A D-optimal design is generated by an iterative search algorithm and seeks to

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
8 of ANSYS, Inc. and its subsidiaries and affiliates.
 Deterministic DoE schemes

minimize the covariance of the parameter estimates for the linear or quadratic approximation model.
This is equivalent to maximizing the determinant , where is the design matrix of model
terms evaluated at all available data points in the design space (Myers and Montgomery 2002).

Remark 1. In optiSLang the number of designs is chosen as 1.5 times the number of Koshal designs. In
LS-OPT, the value for the D-optimality criterion is chosen to be 1.5 times the Koshal design value plus one.

Additionally to the fixed number of designs, optiSLang provides a customizable linear D-optimal design,
which considers a user-defined design number based on a full factorial design for linear regression
functions.

References
Myers, R., and D. C. Montgomery. 2002. Response Surface Methodology. 2nd ed. John Wiley & Sons,
Inc.

2.1.7. Orthogonal Arrays

Available in: LS-OPT, optiSLang

An orthogonal array is a type of fractional factorial design that only considers a selected subset of
the corresponding full factorial samples while still giving complete information about the requested
variable effects.

An orthogonal array is an matrix of experiment, represented by the notation where

is the number of the experiment/design, the number of factors/variables/dimensions and the
number of levels/values. In an orthogonal array, for any pair of columns, all combination of factor
levels occur an equal number of time. This is called the balancing property and it implies orthogon-
ality (Phadke 1989). Table 2.2 (p. 9) provides an orthogonal array example consisting of 7 variables
with 2 levels each.

Table 2.2: orthogonal array which consists of 8 designs in 7-dimensional unit cube

Design
No.
1 0 0 0 0 0 0 0
2 0 0 0 1 1 1 1
3 0 1 1 0 0 1 1
4 0 1 1 1 1 0 0
5 1 0 1 0 1 0 1
6 1 0 1 1 0 1 0
7 1 1 0 0 1 1 0
8 1 1 0 1 0 0 1

Generating the proper orthogonal array suitable for the problem of interest is the main difficulty of
Taguchi method, as a different algorithm may be required depending on the number of samples,
variables and levels. This is the reason why generally, the method uses a pre-identified list of ortho-
gonal arrays.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 9
Design of Experiments

The basic Taguchi method using orthogonal array only considers the main effect of independent
variables. When variable interactions occur, the ability to obtain a good estimation of interaction and
main effects is often lost. In that case, more experiments are needed to estimate all effects. This is
the reason why it is highly recommended to use orthogonal array only in case of non-correlated
variables. Another solution is to use interaction tables. Currently 2-level factors interaction table are
available for and . An interaction table is a square and symmetric matrix where the
dimension is the number of factors. This matrix indicates which column must be left empty in the
orthogonal array to allow the proper analysis of main and interaction effect. For example, if we want
to add an interaction between the factor and the factor using , the column 6 must be
left empty (6 is the value of row 2, column 4 in the interaction table).

Table 2.3: Interaction table for the orthogonal array

Interaction
(1) 3 2 5 4 7 6
(2) 1 6 7 4 5
(3) 7 6 5 4
(4) 1 2 3
(5) 3 2
(6) 1
(7)

To select the appropriate orthogonal array for the factors and levels of interest, the row number of
the matrix of experiments is determined by the matrix degree of freedom. A matrix of experiments
with rows has degrees of freedom. The degree of freedom of an orthogonal array is the sum of
the overall mean degree of freedom (equal to 1), the degree of freedom of every variable and the
degree of freedom of any variable interaction:
(2.1)

where the degree of freedom of a variable is equal to the number of levels minus one. The
degree of freedom of an interaction is equal to the multiplication of both interacting factors degree
of freedom:
(2.2)

As we use a predefined list of orthogonal arrays, it is not always possible to find the perfect orthogonal
array matching the experimental design of interest. When no perfect match can be found, some basic
techniques can be used to select a compatible array while keeping the important orthogonality
property:

• Leave empty columns in a table. A sub-matrix formed by deleting some columns of an orthogonal
array is also an orthogonal array (R. N. Kacker 1991).

• The dummy level technique which assign a factor with levels to a column that has levels with
. The resulting array is still proportionally balanced and hence, orthogonal (Phadke 1989).

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
10 of ANSYS, Inc. and its subsidiaries and affiliates.
 Random and quasi-random DoE schemes

References
Phadke, M. S. 1989. Quality Engineering Using Robust Design. Englewood Cliffs, New Jersey: PTR Prentice-
Hall Inc.

R. N. Kacker, J. J. Filliben, E. S. Lagergren. 1991. "Taguchi"s Orthogonal Arrays Are Classical Designs of
Experiments." Journal of Research of the National Institute of Standards and Technology 96: 577–91.

2.2. Random and quasi-random DoE schemes

2.2.1. Monte Carlo Simulation

Available in: optiSLang

Figure 2.7: Monte Carlo Simulation in a 2-dimensional space.

As alternatives to deterministic design exploration methods, stochastic sampling schemes can be

applied. A very common approach is random sampling, which is called Monte Carlo Simulation (MCS)
(Rubinstein 1981). For design exploration the design variables are assumed to follow a uniform distri-
bution with given lower and upper bounds.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 11
Design of Experiments

The random samples are generated independently in the given design space. If only a small number
of samples is used, often clusters and holes can be observed in the MCS sampling set as shown in
Figure 2.7 (p. 11). More critical is the appearance of undesired correlations between the input variables.
These correlations may have a significant influence on the estimated sensitivity measures.

References
Rubinstein, R. Y. 1981. Simulation and the Monte Carlo Method. New York: John Wiley & Sons.

2.2.2. Latin Hypercube Sampling

Available in: DesignXplorer, LS-OPT, ModelCenter, optiSLang

The Latin Hypercube design (McKay, Beckman, and Conover 1979) is a constrained random experi-
mental design in which, for points, the range of each design variable is subdivided into non-
overlapping intervals on the basis of equal probability. One value from each interval is then selected
at random with respect to the probability density in the interval. The values of the first input
parameter are then paired randomly with the values of the second input parameter. These pairs
are then combined randomly with the values of the third input to form triplets, and so on, until
-tuplets are formed.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
12 of ANSYS, Inc. and its subsidiaries and affiliates.
 Random and quasi-random DoE schemes

Figure 2.8: Latin Hypercube Sampling in a 2-dimensional space.

Latin Hypercube designs are independent of the mathematical model of the approximation and allow
estimation of the main effects of all factors in the design in an unbiased manner. On each level of
every design variable only one point is placed. There are the same number of levels as points, and
the levels are assigned randomly to points. This method ensures that every variable is represented,
no matter if the response is dominated by only a few ones. Another advantage is that the number
of points to be analyzed can be directly defined. Let denote the number of points, and the
number of design variables, each of which is uniformly distributed between 0 and 1. Latin hypercube
sampling (LHS) provides a n-by-k matrix that randomly samples the entire design space broken
down into equal-probability regions:
(2.3)

where are uniform random permutations of the integers 1 through n and are independent
random numbers uniformly distributed between 0 and 1. A common simplified version of LHS has
centered points of equal-probability sub-intervals
(2.4)

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 13
Design of Experiments

LHS can be thought of as a stratified Monte Carlo sampling as shown in Figure 2.8 (p. 13). Latin hy-
percube samples look like random scatter in any bivariate plot, though they are quite regular in each
univariate plot. Often, in order to generate an especially good space filling design, the Latin hypercube
point selection described above is taken as a starting experimental design and then the values in
each column of matrix is permuted so as to optimize some criterion. Several of such criteria are
described in the following.

Figure 2.9: Advanced correlation optimized LHS (left) and space filling LHS using the centered
L2-discrepancy approach (right) in a 2-dimensional space.

References
McKay, M. D., R. J. Beckman, and W. J. Conover. 1979. "A Comparison of Three Methods for Selecting
Values of Input Variables in the Analysis of Output from a Computer Code." Technometrics 21: 239–45.

[Link]. Correlation optimized Latin Hypercube Sampling

Available in: optiSLang

In Advanced Latin Hypercube Sampling, the input distributions and the specified input correlations
are represented very accurately even for a small number of samples. In this approach the spurious
rank order correlations (Iman and Conover 1982) are minimized by a stochastic evolution strategy
(Hungtington and Lyrintzis 1998). For a sufficient performance, the number of samples should be
minimum twice the number of input parameters . Since the pairwise correlations are directly
minimized, the advanced Latin Hypercube Sampling is not limited w.r.t. the number of input para-
meters and can be applied for high dimensional problems.

References
Hungtington, D. E., and C. S. Lyrintzis. 1998. "Improvements to and Limitations of Latin Hypercube
Sampling." Probabilistic Enginerring Mechanics 13: 245–53.

Iman, R. L., and W. J. Conover. 1982. "A Distribution-Free Approach to Inducing Rank Correlation
Among Input Variables." Communications in Statistics - Simulation and Computation 11: 311–34.

[Link]. Space-filling Latin Hypercube Sampling

Available in: DesignXplorer, LS-OPT, optiSLang

One space filling approach is to maximize the minimal distance between any two points,e. g. with
evolution strategies. The strategy would ensure that no two points are too close to each

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
14 of ANSYS, Inc. and its subsidiaries and affiliates.
 Random and quasi-random DoE schemes

other. For a small number of samples , distance designs will generally lie on the exter-
ior of the design space and fill in the interior as becomes larger.

A more robust approach is to minimize the centered L2-discrepancy measure (Hickernell 1998). The
discrepancy is a quantitative measure of non-uniformity of the design points on an experimental
domain. Intuitively, for a uniformly distributed set in the -dimensional cube , we would
expect the same number of points to be in all subsets of having the same volume. Discrepancy
is defined by considering the number of points in the subsets of . Centered L2-discrepancy takes
into account not only the uniformity of the design points over the -dimensional box region ,
but also the uniformity of all the projections of points over lower-dimensional subspaces:

(2.5)

Due to the high computational effort for higher input dimensions, the space-filling Latin Hypercube
sampling is recommended for smaller dimensions with

Remark 1 (Optimized Latin Hypercube sampling). The optiSLang implementations of correlation op-
timized and space-filling Latin Hypercube sampling schemes consider existing start designs in the optim-
ization of the sample permutation.

References
Hickernell, F. J. 1998. "A Generalized Discrepancy and Quadrature Error Bound." Mathematics of
Computation 67: 229–322.

2.2.3. Space-filling designs

Available in: DesignXplorer, LS-OPT, optiSLang

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 15
Design of Experiments

Figure 2.10: Six space-filling designs: 5 points in a 2-dimensional space.

Table 2.4: Description of space-filling algorithms

Algorithm Description
number
0 Random
1 "Central point" Latin Hypercube Sampling (LHS) design with random pairing
2 "Generalized" LHS design with random pairing
3 Given an LHS design, permutes the values in each column of the LHS matrix so as
to optimize the maximin distance criterion taking into account a set of existing
(fixed) design points.
4 Given an LHS design, moves the points within each LHS subinterval preserving the
starting LHS structure, optimizing the maximin distance criterion and taking into
consideration a set of fixed points.
5 Given an arbitrary design (and a set of fixed points), randomly moves the points so
as to optimize the maximin distance criterion using simulated annealing.

In the modeling of an unknown nonlinear relationship, when there is no persuasive parametric regres-
sion model available, and the constraints are uncertain, one might believe that a good experimental
design is a set of points that are uniformly scattered on the experimental domain (design space).
Space-filling designs impose no strong assumptions on the approximation model, and allow a large
number of levels for each variable with a moderate number of experimental points. These designs
are especially useful in conjunction with nonparametric models such as neural networks (section
Section 3.1.9 (p. 40)) and Kriging (section Section 3.1.4 (p. 26)). Space-filling points can also be sub-
mitted as the basis set for constructing an optimal design for a particular model. Some space-filling

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
16 of ANSYS, Inc. and its subsidiaries and affiliates.
 Random and quasi-random DoE schemes

designs are: random Latin Hypercube Sampling (LHS), Orthogonal Arrays, and Orthogonal Latin Hy-
percubes. The key to space-filling experimental designs is in generating "good" random points and
achieving reasonably uniform coverage of sampled volume for a given (user-specified) number of
points. In practice, however, we can only generate finite pseudo-random sequences, which, particularly
in higher dimensions, can lead to a clustering of points, limiting their uniformity. To find a good
space-filling design is a nonlinear programming hard problem, which is difficult to solve exactly. This
problem, however, has a representation, which might be within the reach of currently available tools.
To reduce the search time and still generate good designs, the popular approach is to restrict the
search within a subset of the general space-filling designs. This subset typically has some good "built-
in" properties with respect to the uniformity of a design.

The constrained randomization method termed Latin Hypercube Sampling, has become a popular
strategy to generate points on the "box" (hypercube) design region. The method implies that on each
level of every design variable only one point is placed, and the number of levels is the same as the
number of runs. The levels are assigned to runs either randomly or so as to optimize some criterion,
e.g. so that the minimal distance between any two design points is maximized ("maximin distance"
criterion). Restricting the design in this way tends to produce better Latin hypercubes. However, the
computational cost of obtaining these designs is high. In multidimensional problems, the search for
an optimal Latin hypercube design using traditional deterministic methods may be computationally
prohibitive. This situation motivates the search for alternatives. Probabilistic search techniques, adaptive
simulated annealing and genetic algorithms are attractive heuristics for approximating the solution
to a wide range of optimization problems.

Space-filling designs can be useful for constructing experimental designs for the following purposes:

• The generation of basis points for the D-optimality criterion. This avoids the necessity to create a
very large number of basis points using e.g. the full factorial design for large n. e.g. for n=20 and
3 points per variable, the number of points .

• The generation of design points for all approximation types, but especially for neural networks and
Kriging.

• The augmentation of an existing experimental design. This means that points can be added for
each iteration while maintaining uniformity and equidistance with respect to pre-existing points.

Table 2.4 (p. 16) lists some algorithms to generate space-filling designs.

2.2.4. Sobol Sequences

Available in: optiSLang

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 17
Design of Experiments

Figure 2.11: Sobol sequences in a 2-dimensional space with 63 samples (left)and eight dimensions
with 127 samples (right)

Table 2.5: Recommended number of samples for the Sobol sequences depending on the input
dimension

Dimension Minimum number of samples

63
255
511
1023

The Sobol sequences (Bratley and Fox 1988) is a quasi-random low-discrepancy sampling scheme,
which is a base-2 digital sequence for filling spaces in a highly uniform manner. A full Sobol sequence
consists of samples, whereas can be chosen arbitrary. However, a certain number of samples
is necessary to obtain a space-filling design for higher dimensions, as shown for eight input parameters
and 127 samples in Figure 2.11 (p. 18), where in a specific subset not all quadrants are covered yet.
In Table 2.5 (p. 18) the required number of samples to avoid this phenomena is shown depending
on the input dimension.

References
Bratley, P., and B. L. Fox. 1988. "Algorithm 659 Implementing Sobol"s Quasirandom Sequence Gener-
ator." ACM Transactions on Mathematical Software 14: 88–100.

2.3. Deterministic vs. random DoE schemes

Compared to Latin Hypercube sampling deterministic design schemes have two main disadvantages:
They are limited to a small number of variables due to the rapidly increasing number of required samples
when increasing the model dimension. Further, a reduction of the number of inputs does not improve
the information gained from the samples, since only two or three levels are used in each dimension.
This is illustrated in Figure 2.12 (p. 19). In this figure, a nonlinear function having one major and four

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
18 of ANSYS, Inc. and its subsidiaries and affiliates.
 Glossary

almost unimportant variables is evaluated. Using the LHS, the nonlinearity can be represented very well
in the reduced space. In the case of the full factorial design, which contains three levels in each directions,
again only three positions are left in the reduced space and the dimension reduction does not allow a
better representation of the model response.

Figure 2.12: Dimension reduction for a nonlinear function with five inputs based on Latin
Hypercube Sampling (left) with 100 samples and full factorial design (right) with samples

2.4. Glossary
Bias error. The total error - the difference between the exact and computed response - is composed
of a random and a bias component. The bias component is a systematic deviation between the chosen
model (approximation type) and the exact response of the structure (FEA analysis is usually considered
to be the exact response). Also known as the modeling error. (See also random error).

Center Point. Point at the center value of all factor ranges.

Design of Experiments. See experimental design.

Design matrix. A matrix description of an experiment that is useful for constructing and analyzing ex-
periments.

Design parameter. See design variable.

Design space. A region in the -dimensional space of the design variables ( through ) to which
the design is limited. The design space is specified by upper and lower bounds on the design variables.
Response variables can also be used to bound the design space.

Design variable. An independent design parameter which is allowed to vary in order to change the
design. Symbolized by or (vector containing several design variables).

DOE. Design of Experiments. See experimental design.

D-optimal. The state of an experimental design in which the determinant of the moment matrix
of the least squares formulation is maximized.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 19
Design of Experiments

Effect. How changing the settings of a factor/design variable changes the response. The effect of a
single factor is also called a main effect.

Experimental Design. The selection of designs to enable the construction of a design response surface.
Sometimes referred to as the Point Selection Scheme.

Factors. Process inputs an investigator manipulates to cause a change in the output. Some factors
cannot be controlled by the experimenter but may affect the responses. If their effect is significant,
these uncontrolled factors should be measured and used in the data analysis.

Interaction. Occurs when the effect of one factor on a response depends on the level of another
factor(s).

Iteration. A cycle involving an experimental design, function evaluations of the designs, approximation
and optimization of the approximate problem.

Latin Hypercube Sampling. The use of a constrained random experimental design as a point selection
scheme for response approximation.

Least Squares Approximation. The determination of the coefficients in a mathematical expression so

that it approximates certain experimental results by the minimization of the sum of the squares of the
approximation errors. Used to determine response surfaces as well as calibrating analysis models.

Model. Mathematical relationship which relates changes in a given response to changes in one or more
factors.

Noise. See random error.

Point selection scheme. Same as experimental design.

Process. A series of analysis stages (or steps) designed to produce a result. Multistage process. Example:
metal forming analysis which consists of several stages, e.g. gravity loading, stamping, springback,
trimming, etc.

Random error. The total error - the difference between the exact and computed response - is composed
of a random and a bias component. The random component is, as the name implies, a random deviation
from the nominal value of the exact response, often assumed to be normally distributed around the
nominal value. (See also bias error).

Response. A numerical indicator of the performance of the design.

Saturated design. An experimental design in which the number of points equals the number of unknown
coefficients of the approximation. For a saturated design no test can be made for the lack of fit.

Sampling. In general, Sampling is synonymous with Point Selection or Experimental Design.

Scale factor. A factor which is specified as a divisor of a response in order to normalize the response.

Stochastic. Involving or containing random variables. Involving probability or chance.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
20 of ANSYS, Inc. and its subsidiaries and affiliates.
Chapter 3: Metamodeling techniques
Metamodeling techniques allow the construction of surrogate design models for the purpose of design
exploration such as variable screening, optimization and reliability. The metamodel itself can be used
instead of the numerical or analytical physical model as a fast surrogate model. Before the metamodel
is available, an initial data set, which is often generated with a Design of Experiment (DoE) scheme, has
to be generated with the real physical model with a defined input parameter definition and a suitable
design scheme. Based on the evaluated responses of the physical model, we distinguish between scalar,
signal and field outputs. Once the input and responses data are available, different metmodel types
can be tested and assessed according to the approximation quality. A general framework for an auto-
matic model testing and selection is the Metamodel of Optimal Prognosis in optiSLang. In the following
chapter, first the metamodels for scalar outputs are discussed. Later we present different measures and
procedures for analyzing the model quality, sensitivity evaluation and adaptation strategies.

The following topics are available:

3.1. Regression models for scalar outputs
3.2. Model quality measures
3.3. Model Sensitivity measures
3.4. MOP principle and competition
3.5. Adaptive metamodeling
3.6. Metamodels Glossary

3.1. Regression models for scalar outputs

For scalar outputs, the following meta-models are available within LS-OPT, optiSLang, and DesignXplorer:

• Polynomial regression with linear and quadratic basis functions

• Moving Least Squares (MLS) with linear and quadratic basis functions

• Ordinary Kriging with isotropic and anisotropic kernel functions

• Genetic Aggregation Response Surface (GARS)

• Support Vector Regression (SVR)

• Deep Feed Forward Network (DFFN)

• Deep Infinite Mixture Gaussian Process (DIM-GP)

The polynomial model is the simplest but the fastest available model. It can be applied for small, medium,
and even large number of samples. The approximation of the model e.g. with the MOP solver and
within the exported Functional Mock-Up Unit (FMU) is also faster than all other available models. Moving
Least Squares and Kriging are more time consuming for the model training, especially the anisotropic

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 21
Metamodeling techniques

Kriging, which limits an efficient application up to 1000 samples. You can export these modules as FMU
and are fast in the approximation. MOP filtering is applied to reduced the active variables to the important
ones. Genetic Aggregation Response Surface (GARS) is similar to the anisotropic Kriging limited to
smaller data sets due to exponential increasing training time. Support Vector Regression (SVR) can be
applied even for larger data sets. Both models use the MOP filtering in optiSLang for variable selection.
You currently cannot export FMU for these models. The Deep Feed Forward Network (DFFN) is a deep
learning neural network with automatic feature and variable filtering (Smart Layout). The training of
this model becomes efficient for large data sets with more than 1000 samples, where the training of
other models such as MLS, Kriging, or GARS becomes inefficient. Deep Infinite Mixture Gaussian Process
(DIM-GP) is a further development of the Kriging approach, using a more flexible covariance matrix
description represented by a neural network approximation. The DIM-GP model can be efficiently trained
up to 2000 samples. Currently, variable filtering is not available within the MOP node for the DIM-GP
model. Similar to GARS and SVR, you cannot export the DIM-GP and the DFFN models as FMU.
Table 3.1 (p. 22) describes the properties and a suggested application for each model.

Table 3.1: Available regression models for scalar outputs and recommended application w.r.t.
the number of training data points

Model No. of Variable Recommendation

Samples Filtering
Polynomial arbitrary in MOP Very fast training and fast evaluation
MLS ≤2000 in MOP Fast training and fast evaluation
Isotropic Kriging ≤2000 in MOP Fast training and fast evaluation
Anisotropic ≤500 in MOP Focus on best quality
Kriging
GARS ≤500 in MOP Focus on best quality
SVR ≤5000 in MOP Focus on best quality
DFFN ≥500 Smart layout Works best for large data sets
DIM-GP ≤2000 No Focus on best quality

3.1.1. Polynomial regression

Available in: DesignXplorer, ModelCenter, optiSLang

A commonly used approximation method is polynomial regression, where the model response is
generally approximated by a polynomial basis function of linear or quadratic order with or without
coupling terms.

The model output for a given set of the input parameters can be formulated as the sum of
the approximated value and an error term
(3.1)

where is the polynomial basis,

(3.2)

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
22 of ANSYS, Inc. and its subsidiaries and affiliates.
 Regression models for scalar outputs

and is a vector containing the unknown regression coefficients. These coefficients are generally
estimated from a given set of sampled learning points by assuming independent errors with equal
variance at each point. By using a matrix notation the resulting least squares solution reads
(3.3)

where is a matrix containing the basis polynomials of the learning point samples and is the
vector of learning point values. Further details on the polynomial regression, which is often named
as the classical Response Surface Method can be found in (Montgomery and Runger 2003) and (Myers
and Montgomery 2002).

optiSLang supports linear regression with main terms up to order 2 and linear mixed terms. LS-OPT
allows linear, elliptical (linear and diagonal terms), interaction (linear and off-diagonal terms) and
quadratic functions.

References
Montgomery, D. C., and G. C. Runger. 2003. Applied Statistics and Probability for Engineers. Third. John
Wiley & Sons.

Myers, R., and D. C. Montgomery. 2002. Response Surface Methodology. 2nd ed. John Wiley & Sons,
Inc.

3.1.2. Box-Cox transformation

Available in: DesignXplorer, ModelCenter, optiSLang

Figure 3.1: Box-Cox transformation for different values of

Often the regression model accuracy can be improved by a transformation of the response values.
In (Box and Cox 1964) a family of transformations was introduced:

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 23
Metamodeling techniques

(3.4)

where is geometric mean of the response data and is the transformation parameter. The applic-
ation of the Box-Cox transformation requires the scaling of the response data to pure positive values
e. g. between and . The optimal transformation parameter is obtained by minimizing the re-
gression residuals:

(3.5)

In Figure 3.1 (p. 23) the transformed responses are shown for different values of . The figure indicates
the wide range of possible transformation functions. The Box-Cox transformation was introduced
originally only for linear regression. In optiSlang it is applied for Moving Least Squares and Kriging
as well, by identifying the optimal with a full quadratic regression model.

References
Box, G. E. P., and D. R. Cox. 1964. "An Analysis of Transformations." Journal of the Royal Statistical Society,
Series B 26: 211–52.

3.1.3. Moving Least Squares approximation

Available in: optiSLang

Figure 3.2: Local weighting of support point values of the MLS approximation

In the Moving Least Squares (MLS) approximation (Lancaster and Salkauskas 1981) a local character
of the regression is obtained by introducing position-dependent radial weighting functions. MLS ap-
proximation can be understood as an extension of the polynomial regression.

Similarly the basis function can contain every type of function, but generally only linear and quadratic
terms are used. The approximation function is defined as
(3.6)

with changing ("moving") coefficients in contrast to the constant global coefficients of the
polynomial regression. The final approximation function reads

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
24 of ANSYS, Inc. and its subsidiaries and affiliates.
 Regression models for scalar outputs

(3.7)

where the diagonal matrix contains the weighting function values corresponding to each
learning point. Distance-dependent weighting functions have been introduced. Mostly
the well known Gaussian weighting function is used:

(3.8)

where the influence radius directly influences the width of the weighting function (and the approx-
imation error) and is a numerical constant. A suitable choice of this quantity enables an efficient
smoothing of noisy data. In Figure 3.3 (p. 25) the smoothing effect is shown. In optiSLang the
weighting radius is chosen in the training procedure to get a minimal cross validation error.

Figure 3.3: MLS approximation using a smoothing Gaussian kernel as weighting function (left)
and an interpolating kernel (right) depending on the influence radius

Additional to the classical Gaussian weighting function in Equation 3.8 (p. 25), a regularized weighting
function (Most and Bucher 2005) is available, which enables the direct interpolation of the data points:

(3.9)

where the regularization constant is chosen as . The radius is chosen automatically again
using the cross validation procedure. However, its influence on the approximation function is very
small, as shown in Figure 3.3 (p. 25). The application of this interpolation kernel for noisy data is not
recommended.

References
Lancaster, P., and K. Salkauskas. 1981. "Surface Generated by Moving Least Squares Methods." Math-
ematics of Computation 37: 141–58.

Most, T., and C. Bucher. 2005. "A Moving Least Squares Weighting Function for the Element-Free
Galerkin Method Which Almost Fulfills Essential Boundary Conditions." Structural Engineering and
Mechanics 21: 315–32.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 25
Metamodeling techniques

3.1.4. Kriging
Available in: DesignXplorer, LS-OPT, ModelCenter, optiSLang

Kriging is a method of spatial interpolation composed of two components where the first component
is a polynomial function (representing the "global" trend of the surface) and the second component,
is a realization of a stationary gaussian random process with mean zero and stationary non-negative
covariance (representing the "local" deviations of the surface from the polynomial).

Kriging is named after D. G. Krige (Krige 1951), who applied empirical methods for determining true
ore grade distributions from distributions based on sampled ore grades. In recent years, the Kriging
method has found wider application as a spatial prediction method in engineering design. Detailed
mathematical formulations of Kriging are given by (Forrester, Sobester, and Keane 2008). The basic
postulate of this formulation is:
(3.10)

where is the unknown function of interest, is a known polynomial, the trend model,
and is a stochastic process with mean zero and a given covariance:
(3.11)

With the number of sampling points, is the correlation matrix with the correlation
function between the data points and . is symmetric positive definite with unit diagonal.

The correlation functions are usually defined as exponential type

(3.12)

whereby the exponent is chosen mainly as 1 or 2 which results in an exponential and Gaussian
correlation function, respectively. The unknown correlation lengths can be chosen either as a single
values for all dimensions, which results in an isotropic Kriging approach or individual for all dimensions,
which is called anisotropic Kriging.

The optimal correlation lengths are obtained by a maximum likelihood approach with an optimiz-
ation approach. The final approximation function reads
(3.13)

where is the correlation vector between a prediction point and the support points,
represents the responses at the points and is the trend model. One can choose either a
constant, linear, or quadratic basis function in LS-OPT. The default choice is the constant basis function,
which is also available in optiSLang.

The regression coefficients of the trend model can be derived as

(3.14)

where is a matrix containing the basis polynomials of the support point samples, known from linear
regression. The estimate of the variance of the underlying global model can be derived as

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
26 of ANSYS, Inc. and its subsidiaries and affiliates.
 Regression models for scalar outputs

(3.15)

The presented classical definition of so-called universal Kriging results in an interpolating approximation
function, which respresents the support values directly. In order to get a smoothing effect of the
Kriging, an additional nugget parameter was introduced, which adds to the covariance matrix an in-
dependent noise term. This additional unknown is estimated during the training process within the
maximum likelihood approach (Forrester, Sobester, and Keane 2008). In Figure 3.4 (p. 27) the
smoothing Kriging is shown for a noisy data example.

Figure 3.4: Smoothing Kriging approximation for noisy data

References
Forrester, A., A. Sobester, and A. Keane. 2008. Engineering Design via Surrogate Modelling: A Practical
Guide. John Wiley & Sons.

Krige, D. G. 1951. "A Statistical Approach to Some Basic Mine Valuation Problems on the Witwatersrand."
Journal of the Chemical, Metallurgical and Mining Society of South Africa 52: 119–39.

3.1.5. Radial basis functions

Available in: LS-OPT, ModelCenter, optiSLang

Radial Basis Function (RBF) is an interpolation method that approximates the unknown function
by a linear combinations of radial basis functions centered at given points :

(3.16)

where are the weights of the basis functions and is the basis function value correspond-
ing to center point . In the RBF interpolation approach each of support points is assumed to be
one basis function center point which results in the following formulation:

(3.17)

where the weights of the basis functions can be obtained as:

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 27
Metamodeling techniques

(3.18)

The so-called Gram matrix contains the basis function values of all posible combinations of a pair
of support points and and has to be positive definite. The approximation function of the classical
RBF formulation is a true interpolation which exectly represents the support point values. Typical
basis functions are shown in Figure 3.5 (p. 28). Further details can be found in (Forrester, Sobester,
and Keane 2008). In optiSlang and LS-OPT, Gaussian, multiquadratics, and thin plate spline basis
functions are available.

Figure 3.5: Different Radial Basis Functions formulations depending on the radius

Basis functions
Linear

Cubic

Thin plate
spline
Gaussian

Multiquadratic

Different Radial Basis Functions formulations depending on the radius

In Figure 3.6 (p. 29) the corresponding interpolation functions are shown for the different basis types.
The figure indicates, that all basis type lead to a true interpolation function and all except the linear
type result in a smooth interpolation function.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
28 of ANSYS, Inc. and its subsidiaries and affiliates.
 Regression models for scalar outputs

Figure 3.6: Radial Basis Function interpolation for the different basis types.

Since the classical RBF formulation is a true interpolation, the application for noisy data is not useful.
In (Forrester, Sobester, and Keane 2008) a regularization approach is shown, where the Gram matrix
is modified by adding an additional nugget term to its diagonal
(3.19)

Figure 3.7: Radial Basis Function extension for noisy data.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 29
Metamodeling techniques

In Figure 3.7 (p. 29) the effect of this noise regularization is shown. This smoothing approach is used
in optiSLang together with a polynomial basis model up to fourth order. The optimal scaling parameters
and the nugget parameter are obtained within the training based on the cross validation errors sim-
ilar to the MLS approximation.

If the number of basis function centers are less then the number of supports points, which results in
non-square Gram matrix a least squares estimate of is required, which results in the so-called
Radial Basis Function networks. The center points are typically chosen as a subset of the support
points or arbitrary positions obtain e. g. from cluster analyses.

References
Forrester, A., A. Sobester, and A. Keane. 2008. Engineering Design via Surrogate Modelling: A Practical
Guide. John Wiley & Sons.

3.1.6. Support Vector Regression

Available in: DesignXplorer, LS-OPT, optiSLang

Support Vector Regression (SVR) is a metamodeling technique prescribed for predictably high nonlinear
behavior of the outputs with respect to the inputs. SVR is a non-parametric method based on kernel
functions.

SVR belongs to a general class of Support Vector Method (SVM) type techniques. These are data
classification methods that use hyperplanes to separate data groups. The regression method works
similarly. The main difference is that the hyperplane is used to categorize a subset of the input sample
vectors that are deemed sufficient to represent the output in question. This subset is called the support
vector set. SVR uses non linear kernels.

Given a set of training data that consists of N samples (i=1,...,N) and their corresponding response
values , -SVR regression (Vapnik 2013) aims at finding that has as an upper bound for the
errors on samples while being as smooth as possible. In the case of linear functions, the approximated
function is:
(3.20)

where is a weight vector and is a scalar bias. In order to obtain the "flattest" or most simple
function approximation (to avoid overfitting), the norm of is minimized. Thus the SVR solution is
obtained as the result of the following optimization:

(3.21)

To ensure the feasibility of the optimization, two slack variables and are introduced:

(3.22)

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
30 of ANSYS, Inc. and its subsidiaries and affiliates.
 Regression models for scalar outputs

where is a cost or penalty parameter that ensures a solution with slack variables as close to zero
as possible.

Figure 3.8 (p. 31) provides a simple demonstration.

Figure 3.8: Regression line for a group of sample points with a tolerance of , which is
characterized by slack variables and .

In general, the dual formulation of the SVR optimization is solved:

(3.23)

where and are Lagrange multipliers of the original problem that become the optimization
variables for the dual formulation, is the maximum deviation from the actual response values without
any penalty, are the actual response values, is the cost or penalty parameter and is the number
of training samples. The weight vector eliminated in the dual formulation can be written as:

(3.24)

The support vector regression expression can therefore be written as:

(3.25)

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 31
Metamodeling techniques

where the Lagrange multipliers are obtained by solving the dual problem. The value of can be ob-
tained using the KKT conditions. In the general nonlinear case, the inner product is replaced by a
kernel function :

(3.26)

The optimization is solved using sequential minimal optimization (SMO) (Smola and Schölkopf 2004).
The final SVR expression is:

(3.27)

Note that there are formulations other than the -SVR. For instance, the so-called -SVR (Schölkopf
et al. 2000) controls the number of support vector rather than the errors. There are several possible
settings to train a support vector regression model. We can cite for instance the kernel function, the
parameter value . Among the possible kernels, we can cite polynomial, exponential radial basis
function, Gaussian radial basis function, hyperbolic tangent etc.

The polynomial kernel is given as:

(3.28)

Where is the order of the polynomial.

Remark 1 (Polynomial kernel). In LS-OPT, is internally optimized using 5-fold cross-validation.

The Gaussian RBF kernel is given as:

(3.29)

Remark 2 (Gaussian kernel). Here, is the width parameter or spread of the kernel that is internally op-
timized in LS-OPT using cross-validation. In addition to the kernel parameters, the SVR parameters and
are also internally optimized in LS-OPT during cross-validation.

References
Schölkopf, Bernhard, Alex J Smola, Robert C Williamson, and Peter L Bartlett. 2000. "New Support
Vector Algorithms." Neural Computation 12 (5): 1207–45.

Smola, Alex, and Bernhard Schölkopf. 2004. "A Tutorial on Support Vector Regression." Statistics and
Computing 14 (3): 199–222.

Vapnik, V. 2013. The Nature of Statistical Learning Theory. Springer science & business media.

3.1.7. Sparse Grid Approximation

Available in: DesignXplorer

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
32 of ANSYS, Inc. and its subsidiaries and affiliates.
 Regression models for scalar outputs

The Sparse Grid metamodeling (Garcke 2013) is a hierarchical Sparse Grid interpolation algorithm
based on piecewise multilinear basis functions.

The first ingredient of a Sparse Grid method is a one-dimensional multilevel basis (Figure 3.9 (p. 33)).

Figure 3.9: Piecewise linear hierarchical basis (from the level 0 to the level 3).

The calculation of coefficients values associated to a piecewise linear basis is hierarchical and obtained
by the differences between the values of the objective function and the evaluation of the current
Sparse Grid interpolation (Figure 3.10 (p. 33)).

Figure 3.10: Interpolation with the Hierarchical Basis.

For a multi-dimensional problem, the Sparse Grid metamodeling is based on piecewise multilinear
basis functions that are obtained by a sparse tensor product construction of one-dimensional multi-
level basis.

Let denote a multi-index, where is the problem dimension (number of input

parameters) and corresponds to the level in the i-th direction.

The generation of the Sparse Grid is defined as a sum of grids:

(3.30)

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 33
Metamodeling techniques

In Figure 3.11 (p. 34), an illustration of tensor product in 2D. In Figure 3.12 (p. 35), the combination
of all tensor products in 2D from level 0 to level 3.

Figure 3.11: Tensor product approach to generate the piecewise bilinear basis functions

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
34 of ANSYS, Inc. and its subsidiaries and affiliates.
 Regression models for scalar outputs

Figure 3.12: Tensor product of linear basis functions for a two-dimensional problem

To generate a new Sparse Grid , any Sparse Grid that meets the order relation must
be generated before:

means any tensor product defined in is defined in . For example, in a two-dimensional

problem, the generation of the grid requires several steps: 1. Generation of 2. Generation
of and 3. Generation of 4. Generation of

The calculation of coefficients values associated to a piecewise multilinear basis is similar to the cal-
culation of coefficients of linear basis: the coefficients are obtained by the differences between the
values of the objective function on the new grid and the evaluation (of the same grid) with the current
Sparse Grid interpolation (based on old grids).

You can observe for a higher-dimensional problem that not all input variables carry equal weight. A
regular Sparse Grid refinement can lead to too many support nodes. This is why the Sparse Grid
metamodeling uses a dimension-adaptive algorithm to automatically detect separability and which
dimensions are the more or the less important ones to reduce computational effort for the objectives
functions.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 35
Metamodeling techniques

The hierarchical structure is used to obtain an estimate of the current approximation error. This current
approximation error is used to choose the relevant direction to refine the Sparse Grids. If the approx-
imation error has been found with Sparse Grid , the next iteration consists in the generation of
new Sparse Grids obtained by incrementing of each dimension level of (one by one) as far as
possible: the refinement of can generate two new Sparse Grids and (if the and
already exist).

The Sparse Grid metamodeling stops automatically when the desired accuracy is reached or when
the maximum depth is met in all directions (the maximum depth corresponds to the maximum
number of hierarchical interpolation levels to compute: if the maximum depth is reached in one dir-
ection, the direction is not refined further).

The new generation of the Sparse Grid allows as many linear basis functions as there are points of
discretization (Figure 3.13 (p. 36)).

Figure 3.13: Discretization of piecewise linear basis , and

All Sparse Grids generated by the tensor product contain only one-point that allows refinement more
locally.

In Figure 3.14 (p. 37), Sparse Grid metamodeling is more efficient with a more local refinement process
that uses less sample points, as well as, reaching the requested accuracy faster.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
36 of ANSYS, Inc. and its subsidiaries and affiliates.
 Regression models for scalar outputs

Figure 3.14: Sparse Grid example of 2D response surface.

References
Garcke, Jochen. 2013. "Sparse Grids in a Nutshell." In Sparse Grids and Applications, 57–80. Springer.

3.1.8. Genetic Aggregation Response Surface

Available in: DesignXplorer, optiSLang

Genetic aggregation (Ben Salem and Tomaso 2018) is a metamodeling technique that automatically
selects the best metamodels, including building aggregations of different surrogates.

To select the best metamodel, Genetic Aggregation uses a genetic algorithm that generates populations
of different metamodel solved in parallel. The fitness function of each metamodel is used to determine
which one yields the best approach. It takes into account both the accuracy of the metamodel on
the sample points and the stability of the metamodel (cross-validation).

The Genetic Aggregation response surface (GARS) can be a single metamodel or a combination of
several different metamodels (obtained by a crossover operation during the genetic algorithm).

The Genetic Aggregation response surface can be written as an ensemble using a weighted average
of different metamodels:

(3.31)

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 37
Metamodeling techniques

where: is the prediction of the ensemble. is the prediction of the i-th metamodel. is the
number of metamodels use, is the the weight factor of the i-th metamodel

To estimate the best weight factor values, GARS minimizes the Root Mean Square Error (RMSE) of the
sample points on , and the RMSE of the same points based on the cross-validation ( ).

(3.32)

(3.33)

where:

(3.34)

with:

1. = j-th sample point

2. = the output value for a given sample point

3. = the approximated value of the i-th metamodel built without the j-th sample point

4. ( ) = the number of sample points.

GARS uses the Leave-One-Out and K-Fold cross-validation methods.

Leave-One-Out Method solves as many sub-metamodels as there are sample points. For a given i-th
metamodel, GARS computes sub-metamodels, where each sub-metamodel corresponds to the i-
th metamodel fitted to sample points. The cross-validation error of the j-th sample point is the
error at this point of the sub-metamodel built without the j-th sample point.

K-Fold Method consists in building k sub-metamodels of the i-th metamodel, where each sub-
metamodel corresponds to the i-th metamodel fitted to sample points. The cross-validation
error at the j-th sample point is the error at this point of the sub-metamodel built without the subset
of sample points containing the j-th sample point. To improve the relevance of the k-fold strategy,

the sample points used as validation points of each fold are selected by using the maximum of
minimum distance.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
38 of ANSYS, Inc. and its subsidiaries and affiliates.
 Regression models for scalar outputs

Figure 3.15: Cross-validation error at the second sample point: .

The Leave-One-Out method can be expensive, so GARS uses 10-fold cross-validation by default and
switches to the Leave-One-Out method when the number of sample points is too small to obtain a
relevant 10-fold cross-validation. To reduce the computation cost, the cross-validation is done in
parallel.

If we note , the cross-validation error of i-th metamodel built without the j-th sample point, then:

(3.35)

GARS uses many types of metamodels, including Polynomial Regression (see Section 3.1.1 (p. 22)),
Kriging (see Section 3.1.4 (p. 26)), Support Vector Regression (see Section 3.1.6 (p. 30)), and Moving
Least Squares (see Section 3.1.3 (p. 24)). For each metamodel, there are different settings. For example,
on the Kriging metamodel, you can control the type of kernel (Gaussian, exponential, and so on), the
type of kernel variation (anisotropic or isotropic), and the type of polynomial regression (linear,
quadratic, and so on).

To increase the chance of getting the most effective response surface, GARS generates a population
of metamodels with different types and settings. This population corresponds to the first population
of the genetic algorithm run by GARS. The next populations are obtained by cross-over and mutation
of previous population.

There are two types of cross-over. The first one is the cross-over between two metamodels of the
same type (for example, two Kriging), and the second one is the cross-over between two metamodels
of different types (for example, a Kriging and a polynomial regression).

In the first case, GARS exchanges a part of settings from the first parent to the second parent (for
example, the exchange of kernel type between two Kriging metamodels).

In the second case, GARS creates a new metamodel (an ensemble), which is a combination of the
two parents (for example, the combination of a Kriging and a polynomial regression metamodel).

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 39
Metamodeling techniques

GARS mutates one or several settings of the metamodel (or of the metamodels, in the case of a
combination of metamodel).

To keep a diversity of metamodels, the genetic algorithm removes a part of the metamodel type that
is too present in the population, while it retains other metamodels less present.

In the best case, the population contains similar metamodels in terms of prediction accuracy ( )
when the predicted values are different : this increases the chance of error cancellation on the en-
semble.

References
Ben Salem, M., and L. Tomaso. 2018. "Automatic Selection for General Surrogate Models." Structural
and Multidisciplinary Optimization 58 (2): 719–34.

3.1.9. Neural networks

Available in: LS-OPT, optiSLang

Neural methods are natural extensions and generalizations of regression methods. Neural networks
have been known since the 1940"s, but it took the dramatic improvements in computers to make
them practical (Bishop 1995). Neural networks - just like regression techniques - model relationships
between a set of input variables and an outcome. Neural networks can be thought of as computing
devices consisting of numerical units ("neurons"), whose inputs and outputs are linked according to
specific topologies (see the example in Figure 3.16 (p. 41)). A neural model is defined by its free
parameters - the inter-neuron connection strengths ("weights") and biases. These parameters are
typically "learned" from the training data by using an appropriate optimization algorithm. The training
set consists of pairs of input (design) vectors and associated outputs (responses). The training algorithm
tries to steer network parameters towards minimizing some distance measure, typically the mean
squared error (MSE) of the model computed on the training data.

Several factors determine the predictive accuracy of a neural network approximation and, if not
properly addressed, may adversely affect the solution. For a neural network, as well as for any other
data-derived model, the most critical factor is the quality of training data. In practical cases, we are
limited to a given data set, and the central problem is that of not enough data. The minimal number
of data points required for network training is related to the (unknown) complexity of the underlying
function and the dimensionality of design space. In reality, the more design variables, the more
training samples are required. In the statistical and neural network literature this problem is known
as the "curse of dimensionality". Most forms of neural networks (in particular, feedforward networks)
actually suffer less from the curse of dimensionality than some other methods, as they can concentrate
on a lower-dimensional section of the high-dimensional space. For example, by setting the outgoing
weights from a particular input to zero, a network can entirely ignore that input - see Figure
3.16 (p. 41). Nevertheless, the curse of dimensionality is still a problem, and the performance of a
network can certainly be improved by eliminating unnecessary input variables. It is clear that if the
number of network free parameters is sufficiently large and the training optimization algorithm is
run long enough, it is possible to drive the training MSE error as close as one likes to zero. However,
it is also clear that driving MSE all the way to zero is not a desirable thing to do. For noisy data, this
may indicate over-fitting rather than good modeling. For highly discrepant training data, zero MSE
makes no sense at all. Regularization means that some constraints are applied to the construction of
the neural model with the goal of reducing the "generalization error", that is, the ability to predict
(interpolate) the unobserved response for new data points that are generated by a similar mechanism

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
40 of ANSYS, Inc. and its subsidiaries and affiliates.
 Regression models for scalar outputs

as the observed data. A fundamental problem in modeling noisy and/or incomplete data is to balance
the "tightness" of the constraints with the "goodness of fit" to the observed data. This tradeoff is
called the bias-variance tradeoff in the statistical literature. A multilayer feedforward network and a
radial basis function network are the two most common neural architectures used for approximating
functions. Networks of both types have a distinct layered topology in the sense that their processing
units ("neurons") are divided into several groups ("layers"), the outputs of each layer of neurons being
the inputs to the next layer (Figure 3.16 (p. 41)).

Figure 3.16: Schematic of a neural network with 2 inputs and a hidden layer of 4 neurons with
activation function .

In a feedforward network, each neuron performs a biased weighted sum of their inputs and passes
this value through a transfer (activation) function to produce the output. Activation function of inter-
mediate ("hidden") layers is generally a Sigmoidal function Figure 3.17 (p. 42)A, while network input
and output layers are usually linear (transparent). In theory, such networks can model functions of
almost arbitrary complexity, see (Hornik, Stinchcombe, and White 1990) and (White, Hornik, and
Stinchcombe 1992). All parameters in a feedforward network are usually determined at the same time
as part of a single (non-linear) optimization strategy based on the standard gradient algorithms (the
steepest descent, RPROP, Levenberg-Marquardt, etc.). The gradient information is typically obtained
using a technique called backpropagation, which is known to be computationally effective (Rumelhart,
Hinton, and Williams 1986). For feedforward networks, regularization may be done by controlling the
number of network weights ("model selection"), by imposing penalties on the weights ("ridge regres-
sion") (Hoerl and Kennard 1970), or by various combinations of these strategies (Tikhonov and Arsenin
1977). A radial basis function network has a single hidden layer of radial units, each actually modeling
a response function, peaked at the center, and monotonically varying outwards Figure 3.17 (p. 42)B.
Each unit responds (non-linearly) to the distance of points from its center. The RBF network output
layer is typically linear. Intuitively, it is clear that a weighted sum of the sufficient radial units will always
be enough to model any set of training data (see Figure 3.17 (p. 42)C and Figure 3.17 (p. 42)D). The

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 41
Metamodeling techniques

formal proofs of this property can be found, for example, in (Hartman, Keeler, and Kowalski 1990)
and (Park and Sandberg 1993). An RBF network can be trained extremely quickly, orders of magnitude
faster than a feedforward network. The training process typically takes place in two distinct stages.
First, the centers and deviations of the radial units (i.e. the hidden layer"s weights) must be set; then
the linear output layer is optimized. It is important that deviations are chosen so that RBFs overlap
with some nearby units. Discovering a sub-optimal "spread" parameter typically implies the preliminary
experimental stage. If the RBFs are too spiky, the network will not interpolate between known points
(see Figure 3.17 (p. 42)E). If the RBFs are very broad, the network loses fine detail Figure 3.17 (p. 42)F.
This is actually another manifestation of the over/under-fitting dilemma.

Figure 3.17: Neural Network examples

In the final shape, after training, a multilayer neural network with linear output (Figure 3.16 (p. 41))
can resemble a general linear regression model - a least squares approximation. The major differences
lie in the choice of basis functions and in the algorithms to construct the model (i.e. to adjust model's
free parameters). Techniques to identify the systematical errors in the model and to estimate the
uncertainty of model"s prediction of future observations also become more complex. Unlike polyno-
mial regressors, hidden neurons do not lend themselves to immediate interpretations in terms of input
(design) variables. The next sections discuss various goodness-of-fit assessment approaches applicable
to neural networks. We also discuss how to estimate the variance of the neural model and how to
compute derivatives of a neural network with respect to any of its inputs. Two neural network types,
feedforward and radial basis, are considered.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
42 of ANSYS, Inc. and its subsidiaries and affiliates.
 Regression models for scalar outputs

References
Bishop, C. M. 1995. Neural Networks for Pattern Recognition. Oxford University Press.

Hartman, E. J., J. D. Keeler, and J. M. Kowalski. 1990. "Layered Neural Networks with Gaussian Hidden
Units as Universal Approximations." Neural Computation 2 (2): 210–15.

Hoerl, A. H., and R. W. Kennard. 1970. "Ridge Regression: Biased Estimation for Nonorhtogonal Prob-
lems." Technometrics 12 (3): 55–67.

Hornik, K., M. Stinchcombe, and H White. 1990. "Universal Approximation of an Unknown Mapping
and Its Derivatives Using Multilayer Feedforward Networks." Neural Networks 3: 535–49.

Park, J., and I. W. Sandberg. 1993. "Approximation and Radial Basis Function Networks." Neural Com-
putation 5 (2): 305–16.

Rumelhart, D. E., G. E. Hinton, and R. J Williams. 1986. "Learning Internal Representations by Error
Propagation." Parallel Distributed Processing: Explorations in the Microstructure of Cognition 1: 318–62.

Tikhonov, A. N., and V. Y. Arsenin. 1977. "Solutions of Ill-Posed Problems."

White, H., K. Hornik, and M Stinchcombe. 1992. "Universal Approximation of an Unknown Mapping
and Its Derivatives." Artificial Neural Networks: Approximations and Learning Theory.

[Link]. Feedforward neural networks

Available in: LS-OPT Feedforward (FF) neural networks have a distinct layered topology. Each unit
performs a biased weighted sum of their inputs and passes this value through a transfer (activation)
function to produce the output. The outputs of each layer of neurons are the inputs to the next
layer. In a feedforward network, the activation function of intermediate ("hidden") layers is generally
a sigmoidal function (Figure 3.17 (p. 42))A, network input and output layers being linear. Consider
a FF network with K inputs, one hidden layer with sigmoid units and a linear output unit. For a
given input vector and network weights , the
output of the network is:

(3.36)

Where

The computational graph of Equation 3.36 (p. 43) is shown schematically in Figure 3.16 (p. 41).
The extension to the case of more than one hidden layers can be obtained accordingly. It is
straightforward to show that the derivative of the network Equation 3.36 (p. 43) with respect to
any of its inputs is given by:

(3.37)

Neural networks have been mathematically shown to be universal approximators of continuous

functions and their derivatives (on compact sets) (Hornik, Stinchcombe, and White 1990). In other
words, when a network (White, Hornik, and Stinchcombe 1992) converges towards the underlying
function, all the derivatives of the network converge towards the derivatives of this function.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 43
Metamodeling techniques

Standard non-linear optimization techniques including a variety of gradient algorithms (the steepest
descent, RPROP, Levenberg-Marquardt, etc.) are applied to adjust FF network's weights and biases.
For neural networks, the gradients are easily obtained using a chain rule technique called 'back-
propagation' (Rumelhart, Hinton, and Williams 1986). The second-order Levenberg-Marquardt al-
gorithm appears to be the fastest method for training moderate-sized FF neural networks (up to
several hundred adjustable weights) (Bishop 1995). However, when training larger networks, the
first-order RPROP algorithm becomes preferable for computational reasons (Riedmiller and Braun
1993).

Regularization: For FF networks, regularization may be done by controlling the number of network
weights ('model selection'), by imposing penalties on the weights ('ridge regression'), or by various
combinations of these strategies ((Hoerl and Kennard 1970), (Tikhonov and Arsenin 1977)). Model
selection requires choosing the number of hidden nodes and, sometimes, the number of network
hidden layers. Most straightforward is to search for an 'optimal' network architecture that minimizes
, or (see further details in Section 3.2.2 (p. 58) and Section 3.2.3 (p. 58)).
Often, it is feasible to loop over 1, 2,... hidden nodes and finally select the network with the smallest
GCV error. In any event, in order for the GCV measure to be applicable, the number of training
points should not be too small compared to the required network size M.

Over-fitting: To prevent over-fitting, it is always desirable to find neural solutions with the smallest
number of parameters. In practice, however, networks with a very parsimonious number of weights
are often hard to train. The addition of extra parameters (i. e. degrees of freedom) can aid conver-
gence and decrease the chance of becoming stuck in local minima or on plateaus (Lawrence, Giles,
and Tsoi 1996). Weight decay regularization involves modifying the performance function , which
is normally chosen to be the mean sum of squares of the network errors on the training set
(Equation 3.51 (p. 55)). When minimizing MSE (Equation 3.51 (p. 55)) the weight estimates tend to
be exaggerated. We can impose a penalty for this tendency by adding a term that consists of the
sum of squares of the network weights (see also Equation 3.51 (p. 55)):
(3.38)

where

(3.39)

where is the number of weights and the number of points in the training set. Notice that
network biases are usually excluded from the penalty term . Using the modified performance
function (Equation 3.38 (p. 44)) will cause the network to have smaller weights, and this will force
the network response to be smoother and less likely to overfit. This eliminates the guesswork required
in determining the optimum network size. Unfortunately, finding the optimal value for and is
not a trivial task. If we make is too small, we may get over-fitting. If is too large, the
network will not adequately fit the training data. A rule of thumb is that a little regularization usually
helps (Sjöberg and Ljung 1992). It is important that weight decay regularization does not require
that a validation subset be separated out of the training data. It uses all of the data. This advantage
is especially noticeable in small sample size situations. Another nice property of weight decay reg-
ularization is that it can lend numerical robustness to the Levenberg-Marquardt algorithm. The L-
M approximation to the Hessian of Equation 3.38 (p. 44) is moved further away from singularity
due to a positive addend to its diagonal:
(3.40)

where

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
44 of ANSYS, Inc. and its subsidiaries and affiliates.
 Regression models for scalar outputs

(3.41)

In (Bishop 1995), (MacKay 1992), (Foresee and Hagan 1997) and (Moody 1992) the Bayesian ('evidence
framework' or 'type II maximum likelihood') approach to regularization is discussed. The Bayesian
re-estimation algorithm is formulated as follows. At first, we choose the initial values for and .
Then, a neural network is trained using a standard non-linear optimization algorithm to minimize
the error function (Equation 3.38 (p. 44)). After training, i. e. in the minimum of Equation 3.38 (p. 44),
the values for and are re-estimated, and training restarts with the new performance function.
Regularization hyper-parameters are computed in a sequence of 3 steps:

(3.42)

where are (positive) eigenvalues of matrix in Equation 3.40 (p. 44), is the estimate
of the effective number of parameters of a neural network,
(3.43)

It should be noted that the algorithm (Equation 3.42 (p. 45)) relies on numerous simplifications
and assumptions, which hold only approximately in typical real-world problems (Cohn 1996). In
the Bayesian formalism a trained network is described in terms of the posterior probability distribu-
tion of weight values. The method typically assumes a simple Gaussian prior distribution of weights
governed by an inverse variance hyper-parameter . If we present a new input vector
to such a network, then the distribution of weights gives rise to a distribution of network outputs.
There will be also an addend to the output distribution arising from the assumed
Gaussian noise on the output variables:
(3.44)

With these assumptions, the negative log likelihood of network weights given training points
x(1), ... , x(N) is proportional to MSE (Equation 3.51 (p. 55)), i. e., the maximum likelihood estimate
for is that which minimizes (Equation 3.51 (p. 55)) or, equivalently, . In order for Bayes estimates
of and to do a good job of minimizing the generalization in practice, it is usually necessary
that the priors on which they are based are realistic. The Bayesian formalism also allows us to cal-
culate error bars on the network outputs, instead of just providing a single 'best guess' output .
Given an unbiased model, minimization of the performance function (Equation 3.51 (p. 55)) amounts
to minimizing the variance of the model. The estimate for output variance of the network at
a particular point is given by:
(3.45)

Equation 3.45 (p. 45) is based on a second-order Taylor series expansion of Equation 3.38 (p. 44)
around its minimum and assumes that is locally linear. Neural networks have a natural
variability because of the following reasons:

• Local behavior of the neural network training algorithms

• Uncertainty (noise) in the training data.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 45
Metamodeling techniques

The neural network training error function usually has multiple local and global minima. With dif-
ferent initial weights, the training algorithm typically ends up in different (but usually almost equally
good/bad) local minima. The larger the amount of noise in the data, the larger is the difference
between these NN solutions.

References
Bishop, C. M. 1995. Neural Networks for Pattern Recognition. Oxford University Press.

Cohn, David A. 1996. "Neural Network Exploration Using Optimal Experiment Design." Neural Networks
9 (6): 1071–83.

Foresee, F Dan, and Martin T Hagan. 1997. "Gauss-Newton Approximation to Bayesian Learning."
In Proceedings of International Conference on Neural Networks (ICNN'97), 3:1930–35. IEEE.

Hoerl, A. H., and R. W. Kennard. 1970. "Ridge Regression: Biased Estimation for Nonorhtogonal
Problems." Technometrics 12 (3): 55–67.

Hornik, K., M. Stinchcombe, and H White. 1990. "Universal Approximation of an Unknown Mapping
and Its Derivatives Using Multilayer Feedforward Networks." Neural Networks 3: 535–49.

Lawrence, S. C., Lee Giles, and Ah Chung Tsoi. 1996. "What Size Neural Network Gives Optimal
Generalization? Convergence Properties of Backpropogation." Technical Report UMIACS-TR-96-22
and CS-TR-3617. University of Maryland.

MacKay, David JC. 1992. "Bayesian Interpolation." Neural Computation 4 (3): 415–47.

Moody, J. E. 1992. "The Effective Number of Parameters: An Analysis of Generalization and Regular-
ization in Nonlinear Learning Systems." Neural Computation 4 (3): 415–47.

Riedmiller, M., and H. Braun. 1993. "A Direct Adaptive Method for Faster Backpropagation Learning:
The RPROP Algorithm." Proceedings of the IEEE International Conference on Neural Networks (ICNN)
1: 586–91.

Rumelhart, D. E., G. E. Hinton, and R. J Williams. 1986. "Learning Internal Representations by Error
Propagation." Parallel Distributed Processing: Explorations in the Microstructure of Cognition 1: 318–62.

Sjöberg, J., and L. Ljung. 1992. "Overtraining, Regularization, and Searching for Minimum in Neural
Networks." Preprints of the 4th IFAC Int. Symp. On Adaptive Systems in Control and Signal Processing,
669.

Tikhonov, A. N., and V. Y. Arsenin. 1977. "Solutions of Ill-Posed Problems."

White, H., K. Hornik, and M Stinchcombe. 1992. "Universal Approximation of an Unknown Mapping
and Its Derivatives." Artificial Neural Networks: Approximations and Learning Theory.

[Link]. Radial basis function networks

Available in: LS-OPT

A radial basis function (Equation 3.16 (p. 27)) can be interpreted as an artificial radial basis function
called radial basis function neural network. A radial basis function neural network has a distinct 3-
layer topology. The input layer is linear (transparent). The hidden layer consists of non-linear radial

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
46 of ANSYS, Inc. and its subsidiaries and affiliates.
 Regression models for scalar outputs

units, each responding to only a local region of input space. The output layer performs a biased
weighted sum of these units and creates an approximation of the input-output mapping over the
entire space.

The activation of radial basis function is determined by the Euclidean distance

between the input vector and RBF center in -dimensional space.

The Gaussian basis function is a localized function (peaked at the center and descending outwards)
with the property that as . Parameter controls the smoothness properties of the RBF
unit.

For a given input vector x the output of RBF network with K inputs and a hidden layer with H basis
function units is given by:

(3.46)

where

(3.47)

Notice that hidden layer parameters represent the center of radial unit, while
corresponds to its deviation. Parameters and are the output layer"s bias and
weights, respectively. A linear super-position of localized functions as in Equation 3.38 (p. 44) is
capable of universal approximation. The formal proofs of this property can be found, for example,
in (Hartman, Keeler, and Kowalski 1990) and (Park and Sandberg 1993). It is straightforward to show
that the derivative of the network Equation 3.38 (p. 44) with respect to any of its inputs is given
by:

(3.48)

where denotes the first derivative of the transfer function .

Theory tells us that when a network (Equation 3.46 (p. 47)) converges towards the underlying
function, all the derivatives of the network converge towards the derivatives of this function. A key
aspect of RBF networks, as distinct from feedforward neural networks, is that they can be interpreted
in a way which allows the hidden layer parameters (i.e. the parameters governing the radial functions)
to be determined by semi-empirical, unsupervised training techniques. Accordingly, although a
radial basis function network may require more hidden nodes than a comparable feedforward
network, RBF networks can be trained extremely quickly, orders of magnitude faster than FF networks.

For RBF networks, the training process generally takes place in two distinct stages. First, the centers
and deviations of the radial units (i.e. hidden layer parameters and ) must
be set. All the basis functions are then kept fixed, while the linear output layer (i.e., ) is
optimized in the second phase of training. In contrast, all of the parameters in a FF network are
usually determined at the same time as part of a single training (optimization) strategy. Techniques
for selecting and are discussed at length in following paragraphs. Here

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 47
Metamodeling techniques

we shall assume that the RBF parameters have already been chosen, and we focus on the problem
of optimizing the output layer weights.

Mathematically, the goal of output layer optimization is to minimize a performance function, which
is normally chosen to be the mean sum of squares of the network errors on the training set
(Equation 3.51 (p. 55)). If the hidden layer"s parameters in Equation 3.37 (p. 43)
are kept fixed, then the performance function (Equation 3.51 (p. 55)) is a quadratic function of the
output layer" parameters and its minimum can be found in terms of the solution of a
set of linear equations (e.g., using singular value decomposition). The possibility of avoiding the
need for time-consuming and costly non-linear optimization during training is one of the major
advantages of RBF networks over FF networks. However, when the number of optimized parameters
( , in our case) is small enough, non-linear optimization (Levenberg-Marquardt, etc.) may also
be cost-effective.

It is clear that the ultimate goal of RBF neural network training is to find a smooth mapping which
captures the underlying systematic aspects of the data without fitting the noise. However, for noisy
data, the exact RBF network, which passes exactly through every training data point, is typically a
highly oscillatory function. There are a number of ways to address this problem. By analogy with
FF network training, one can add to Equation 3.51 (p. 55) a regularization term that consists of the
mean of the sum of squares of the optimized weights. In conventional curve fitting this form of
regularization is called ridge regression. The sub-optimal value for hyperparameters and in
Equation 3.38 (p. 44) can be found by applying Bayesian re-estimation formulae (Equation
3.40 (p. 44)) - (Equation 3.42 (p. 45)). It is also feasible to iterate over several trial values of and
.

For RBF networks, however, the most effective regularization methods are probably those pertaining
to selecting radial centers and deviations in the first phase of RBF training. The commonly held
view is that RBF centers and deviations should be chosen so as to form a representation of the
probability density of the input data. The classical approach is to set RBF centers equal to all the
input vectors from the training dataset. The width parameters are typically chosen - rather arbit-
rarily - to be some multiple of the average spacing between the RBF centers (e.g. to be roughly
twice the average distance). This ensures that the RBF"s overlap to some degree and hence give a
relatively smooth representation of data.

To simplify matters, the same value of the width parameter for all RBF units is usually considered.
Sometimes, instead of using just one value for all RBF"s, each RBF unit"s deviation is individually
set to the distance to its nearest neighbors. Hence, deviations become smaller in densely
populated areas of space, preserving fine detail, and are higher in sparse areas of space, interpolating
between points where necessary. Again the choice of is somewhat arbitrary. RBF networks with
individual radial deviations can be particularly useful in situations where data tend to cluster in
only a small subregion of the design space (for example, around the optimum of the underlying
system which RSM is searching for) and are sparse elsewhere.

One must take into consideration that after the first phase of RBF training is over, there"s no way
to compensate for large inaccuracies in radial deviations by, say, adding a regularization term
to the performance function. If the basis functions are too spiky, the network will not interpolate
between known points, and thus, will lose the ability to generalize. If the Gaussians are very broad,
the network is likely to lose fine detail. The popular approach to find a sub-optimal spread parameter
is to loop over several trial values of and , and finally select the RBF network with the smallest
GCV (FPE, CV-k) error. This is somewhat analogous to searching for an optimal number of hidden
nodes of a feedforward neural network.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
48 of ANSYS, Inc. and its subsidiaries and affiliates.
 Regression models for scalar outputs

In order to eliminate all the guesswork required in determining RBF deviations, it might seem nat-
ural to treat ( , to be precise) in Equation 3.46 (p. 47) as adjustable para-
meters, which are optimized in the second phase of training along with the output layer"s weights
and bias. Practical applications of this approach, however, are rare. The reason may be that it requires
the use of a non-linear optimization method in combination with a sophisticated regularization
scheme specially designed so as to guarantee that the radial functions will remain localized.

It should be noted that RBF networks may have certain difficulties if the number of RBF units ( )
is large. This is often the case in multidimensional problems. The difficulty arises because for a large
number of RBF"s, a large number of training samples are required in order to ensure that the
neural network parameters are properly determined. A large number of RBF units also increase the
computation time spent on optimization of the network output layer and, consequently, the RBF
architecture loses its main (if not the only one) advantage over FF networks - fast training.

Radial basis function networks actually suffer more from the curse of dimensionality than feedforward
neural networks. To explain this statement, consider the effect of adding an extra, perfectly spurious
input variable to a network. A feedforward network can learn to set the outgoing weights of the
spurious input to zero, thus ignoring it. An RBF network has no such luxury: data in the relevant
lower-dimensional space get "smeared" out through the irrelevant dimension, requiring larger
numbers of units to encompass the irrelevant variability.

In principle, the number of RBF"s ( ) need not equal the number of training samples ( ), and RBF
units are not constrained to be centered on the training data points. In fact, when data contain
redundant information, we do not need all data points in learning. One simple procedure for selecting
RBF centers is to set them equal to a random subset of the input vectors from the training set.
Since they are randomly selected, they will "represent" the distribution of the (redundant) training
data in a statistical sense. Of course, and should not be too small in this case.

It is clear, however, that the optimal choice of RBF centers based on the input data alone need not
be optimal for representing the input-output mapping as reflected in the observed data. In order
to overcome these limitations, the selection procedure should take into account the output values,
or at least, approximate estimates (assumptions) of the global behavior of the underlying system.
The common neural term for such techniques involving output values is "active learning". In the
context of active learning, RBF networks can be thought of as DOE metamodels analogous to
polynomials, (MacKay 1992) and (Cohn 1996). Given a candidate list of points, an active learner is
searching for the "best" points in order to position RBF centers. Popular in neural applications is to
treat RBF active learning as "pruning" technique intended for identifying critical data and discarding
redundant points, or more accurately, not selecting some training points as RBF centers. RBF active
learning methods are being successfully applied to approximate huge datasets that come from
natural stochastic processes. It is questionable, however, whether active learning can be useful for
non-redundant datasets, specifically for RSM design sets generated by performing DOE analysis
based on low-order polynomial metamodels.

To briefly summarize, parameters governing radial units (radial centers and deviations) play a key
role in generalization performance of a RBF model. The appropriate selection of RBF centers implies
that we choose a minimal number of training data points that carry enough information to build
an adequate input-output representation of the underlying function. Unfortunately, this is easier
said than done. Indeed, there is a general agreement that selecting RBF centers and deviations is
more Art than Science.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 49
Metamodeling techniques

References
Cohn, David A. 1996. "Neural Network Exploration Using Optimal Experiment Design." Neural Networks
9 (6): 1071–83.

Hartman, E. J., J. D. Keeler, and J. M. Kowalski. 1990. "Layered Neural Networks with Gaussian Hidden
Units as Universal Approximations." Neural Computation 2 (2): 210–15.

MacKay, David JC. 1992. "Bayesian Interpolation." Neural Computation 4 (3): 415–47.

Park, J., and I. W. Sandberg. 1993. "Approximation and Radial Basis Function Networks." Neural
Computation 5 (2): 305–16.

[Link]. Deep Feed Forward Network

Available in: optiSLang

Deep Feed Forward Network (DFFN) implements a metamodel based on Artificial Neural Network
paradigms. It uses the Keras library with TensorFlow as backend to build and train the network.

Compared to traditional feedforward networks proposed in section Section [Link] (p. 43), the Deep
Feed Forward Network (DFFN) contains multiple hidden layers between the input and output layers,
see Figure 3.18 (p. 50), which makes it more capable of learning intricate patterns and representations
in data.

Figure 3.18: Schematic view of a Deep Feed Forward Network with 2 inputs and 3 hidden
layers of 4 neurons with activation function f.

The identification of the network structure and respective hyper parameters is an important pre-
requisite to achieve high model accuracy during model training. Therefore, the DFFN comes with
different network layout options. Depending on the layout chosen, the construction of the neural
network model involves a different strategy for hyperparameter search as well as model training,
by means of identifying the optimal network weights.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
50 of ANSYS, Inc. and its subsidiaries and affiliates.
 Regression models for scalar outputs

[Link].1. Hyperparameter search

Hyperparameter search is the process of finding the optimal topology and configuration for a
neural network to be trained. An optimal configuration of the hyperparameter is crucial for the
resulting approximation quality. With DFFN, two strategies are provided to support an automatic
identification of the optimal hyperparameter configuration, namely smart layout and optimization
layout. This comprises the identification of the following Hyperparameters:

Learning rate

Controls how much to change the model in response to the estimated error each time the
model weights are updated.

Number of Hidden Layers and Units

Determines the architecture of the neural network.

Activation Functions

Functions that determine the output of a neural network node.

Regularization coefficient

A parameter that controls the amount of regularization applied to the model to prevent
over-fitting. In the L1 regularization technique used in the Deep Feed Forward Network, the
absolute value of weights is added to the loss function:

(3.49)

where is the original loss, is the regularization coefficient, and are the weights. It can
lead to sparse models where individual weights will be set to zero.

In the manual layout, all Hyperparameters are to be configured by the user, requiring expert
knowledge on the training of artificial neural networks.

[Link].2. Model training

The training of the neural network is carried out in a cross-validation setup. The cross-validation
is a technique used in all network layouts to evaluate the performance of deep feed forward
networks, by partitioning the data into subsets and ensuring a high level of robustness and
generalizability of the model. The trained models are then used to obtain a set of predictions.
The training steps are described below:

Data split

The dataset is divided into subsets or "folds". Common choices for are between
to .

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 51
Metamodeling techniques

Training and validation

The model is trained times. Each time, one of the subsets is used as the validation
set, and the remaining subsets are used as the training set.

Predictions averaging

The predictions are calculated for each of the iterations. The final prediction is the average
of these values.

This process helps in detecting over-fitting and ensures that the model performs well on unseen
data.

Figure 3.19: Schematic view on cross-validation procedure used for the training and
evaluation of the DFFN model.

[Link].3. Network layout

[Link].3.1. Smart layout

Smart layout refers to a smart process recommended for first-time user and for non machine
learning experts. It automatically adjusts the hyperparameters to ensure ease of use and quality
of prediction. Only parameters related to model architecture (maximum hidden layers and
maximum number of neurons per layer) are user-defined. In this process, an input sensitivity
and feature selection are performed at the end of the training. It is used to determine the most
influential neurons and to reduce the input space accordingly. This reduction takes into account
its impact on the CoP (section Section 3.2.5 (p. 61)), given a respective tolerance value. A post-
training pruning is applied to remove individual neurons based on their importance. This al-
gorithmic step consists of identifying and removing the neurons, which contribute the least to
the model performance. The reduced network model is more compact, consumes less memory,

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
52 of ANSYS, Inc. and its subsidiaries and affiliates.
 Regression models for scalar outputs

requires less computational power for approximation and most relevant is expected to gener-
alize better on unseen data.

[Link].3.2. Optimization layout

The optimization option gives experienced users greater flexibility in selecting the hyper-
parameters they wish to integrate into the search hyperparameter. A large number of parameters
can be set, impacting on model architecture, the training process, etc..

[Link].3.3. Manual layout

Compared with the two previous options, the manual layout does not incorporate a hyperpara-
meter search. The user can manually select parameter values to define network architecture
and training parameters. This layout is recommended for users with a good level of expertise
in deep learning.

Remark 1. Using Deep Feed Forward Network for creating a MOP is beneficial if you have a lot of
input designs, as its required computational time scales better than that of other surrogate models.

3.1.10. Deep Infinite Mixture Gaussian Process (DIM-GP)

Available in: optiSLang

Deep Infinite Mixture of Gaussian Processes (DIM-GP) (Cremanns et al. 2021) is a probabilistic machine
learning algorithm consisting of a unique combination of neural networks (Goodfellow, Bengio, and
Courville 2016) and Gaussian processes (Williams and Rasmussen 2006).

Figure 3.20: Schematic representation of DIM-GP

This combines many of the advantages of both methods and thus negates some of their disadvantages.
In most machine learning models, the so-called free parameters are fixed after training and do not
change. With DIM-GP, these parameters are dependent on the point of prediction. Thus, the model
behavior of DIM-GP can change within a design space as if several models had been trained on dif-
ferent areas of the design space. This makes DIM-GP particularly flexible and allows it to better rep-
resent non-continuous output variables or output variables where the physical behavior changes.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 53
Metamodeling techniques

This flexibility is achieved by using a neural network to predict the free parameters of a Gaussian
process depending on input parameter values. i.e. when training DIM-GP a neural network is trained
but in the end a Gaussian process is used for the prediction. Another neural network is used to learn
the influence of noise in the data, since this is also a free parameter in the Gaussian process. Here,
the noise is also predicted by the neural network depending on the values of the input parameters
and thus is not fixed after training. Further properties of DIM-GP are:

• Applicability of Gaussian processes even with many samples by mini-batch training. i.e. not all
samples are loaded into the working memory at the same time, but processed batchwise. This is
possible because neural networks are trained in DIM-GP. Indirectly this means that during mini-
batch training new Gaussian processes are formed again and again, and this corresponds to a
training of an ensemble of "infinitely" many Gaussian processes.

• No setting of hyperparameters necessary. In the Gaussian process, there are basically no hyperpara-
meters except the choice of the covariance function. Depending on the relationship between the
output and input parameters, certain covariance functions may be more appropriate than others.
In addition, in certain cases a so-called kernel engineering is advantageous where more than one
covariance function is used at the same time. e.g., in order to be able to represent different super-
imposed physical effects better. DIM-GP takes over the choice of the covariance function independ-
ently or performs an automatic kernel engineering, which is even dependent on the values of the
input parameters (non-stationary). As a result, DIM-GP typically does not require user settings that
could have a large impact on the model. The neural networks in DIM-GP are fixed and would not
need to be changed.

• Automatic detection of noise and outliers. Since Gaussian processes have the ability to handle noisy
data, DIM-GP also has this feature. DIM-GP does not need to repeat experiments but learns this
automatically. This is an option that the user can turn on or off.

• Applicable to regression and classification tasks.

• Applicable to scalar output quantities as well as signals.

References
Cremanns, Kevin, Dirk Roos, Stefan Reh, and Sebastian Münstermann. 2021. "Probabilistic Machine
Learning for Pattern Recognition and Design Exploration." Fachgruppe für Materialwissenschaft und
Werkstofftechnik.

Goodfellow, Ian, Yoshua Bengio, and Aaron Courville. 2016. Deep Learning. MIT press.

Williams, Christopher KI, and Carl Edward Rasmussen. 2006. Gaussian Processes for Machine Learning.
Vol. 2. 3. MIT press Cambridge, MA.

3.2. Model quality measures

Response surfaces are useful for interactive trade-off studies. For the trade-off surface to be useful, its
capability of predicting accurate response must be known. Error analysis is therefore an important aspect
of using response surfaces in design. Inaccuracy is not always caused by random error (noise) only, but
modeling error (sometimes called bias error), especially in a large subregion or where there is strong
non-linearity present, could play a very significant role. There are several error measures available to
determine the accuracy of a response surface.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
54 of ANSYS, Inc. and its subsidiaries and affiliates.
 Model quality measures

3.2.1. Residual analysis

[Link]. Residual sum of squares

Available in: DesignXplorer, LS-OPT

For the predicted response and the actual response , this error is expressed as:

(3.50)

If applied only to the regression points, this error measure is not very meaningful unless the design
space is oversampled e.g., = 0 if the number of points equals the number of basis functions L
in the approximation.

[Link]. Mean squared error (MSE)

Available in: LS-OPT

The Mean squared error is given by:

(3.51)

[Link]. Root mean squared error (RMSE)

Available in: DesignXplorer, LS-OPT, optiSLang

The residual sum-of-squares is sometimes used in its square root form, , and called the RMS
error (so called RMSE or Root Mean Squared Error)

(3.52)

[Link]. Relative root mean squared error

Available in: DesignXplorer

This is the square root of the average squared of the residuals scaled by the actual output values
at the sample points for regression methods. The best value is 0. In general, the closer the value is
to 0, the better quality of the response surface.

However, in some situations, you can have a larger value and still have a good response surface.
For example, this can be true when some of the output values are close to zero (the observed value
= 1e-10 while the predicted value = 1.0e-8).

(3.53)

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 55
Metamodeling techniques

[Link]. Maximum residual

Available in: LS-OPT, optiSLang

This is the maximum residual considered over all the design points and is given by:

(3.54)

[Link]. Maximum Relative Residual

Available in: DesignXplorer

The maximum relative residual is the maximum distance out of all generated points from the cal-
culated response surface to each generated point. In general, the closer the value is to 0, the better
quality of the response surface. However, in some situations, you can have a larger value and still
have a good response surface. For example, this can be true when the mean of the output values
is close to zero.

(3.55)

[Link]. Relative Maximum Absolute Error

Available in: DesignXplorer

The relative maximum absolute error is the absolute maximum residual value relative to the
standard deviation of the actual output data, modified by the number of samples. The best value
is 0. In general, the closer the value is to 0, the better quality of the response surface.

The relative maximum absolute error is mathematically represented as:

(3.56)

[Link]. Relative Average Absolute Error

Available in: DesignXplorer

The relative average absolute error is the average of the residuals relative to the standard deviation
of the actual outputs. This value is useful when the number of samples is low ( ). The best value
is 0. In general, the closer the value is to 0, the better quality of the response surface.

The relative average absolute error and the relative maximum absolute error correspond to the
maximum error and average absolute error scaled by the standard deviation. For example, the rel-
ative root mean square error becomes negligible if both of these values are small.

The relative average absolute error is mathematically represented as:

(3.57)

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
56 of ANSYS, Inc. and its subsidiaries and affiliates.
 Model quality measures

[Link]. PRESS residuals for polynomial regression

Available in: LS-OPT

The prediction sum of squares residual (PRESS) uses each possible subset of -1 points as a regres-
sion data set to fit the regression model, and the remaining point in turn is used to form a prediction
set (Myers and Montgomery 2002). PRESS is defined as the sum of all prediction errors, but can be
computed from a single regression analysis of all points.

(3.58)

where are the diagonal terms of

(3.59)

where is a matrix containing the basis polynomials of the learning point samples and is the
"hat" matrix, the matrix that maps the observed responses to the fitted responses, i.e.
(3.60)

The PRESS residual can also be written in its square root form

(3.61)

For a saturated design, equals the unit matrix so that the PRESS indicator becomes undefined.

References
Myers, R., and D. C. Montgomery. 2002. Response Surface Methodology. 2nd ed. John Wiley & Sons,
Inc.

[Link]. R2 for prediction of polynomial regression

Available in: LS-OPT

For the purpose of accuracy the indicator has been devised (Myers and
Montgomery 2002).
(3.62)

where

(3.63)

where represents the ability of the model to detect the variability in predicting new re-
sponses.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 57
Metamodeling techniques

References
Myers, R., and D. C. Montgomery. 2002. Response Surface Methodology. 2nd ed. John Wiley & Sons,
Inc.

3.2.2. Generalized cross-validation and Akaike"s final prediction error

Available in: LS-OPT The generalized cross-validation (GCV) (Wahba 1990) and Akaike"s final prediction
error (FPE) (Akaike 1998) provide computationally feasible means of estimating the appropriateness
of the model. GCV and FPE estimates combine the training MSE with a measure of the model com-
plexity:
(3.64)

(3.65)
(3.66)
(3.67)

where is the (effective) number of model parameters. In theory, GCV estimates should be related
to . As a very rough approximation to , we can assume that all of the network free parameters are
well determined so that = M, where M is the total number of network weights and biases. This is
what we would expect to be the case for large so that . Note that GCV is undefined when
is equal to the number of training points (N). In theory, GCV and FPE estimates should be related
to the effective number of model"s parameters . Techniques to assess for neural networks will be
discussed later. GCV and FPE measures are asymptotically equivalent for large N.

References
Akaike, Hirotugu. 1998. "Statistical Predictor Identification." Selected Papers of Hirotugu Akaike, 137–51.

Wahba, Grace. 1990. Spline Models for Observational Data. SIAM.

3.2.3. k-fold cross-validation

Available in: DesignXplorer, LS-OPT, optiSLang

The -fold cross-validation (denoted here as CV- ) provides computationally feasible means of estim-
ating the appropriateness of the model. In -fold cross-validation the training dataset is divided into
randomly selected disjoint subsets of roughly equal size . The model is trained and tested
times. Each time it is trained on all data except for points from subset and then tested on -th
subset. Formally, let be the prediction of such a model for the points from
subset . Then the CV-k estimates of accuracy:

(3.68)

(3.69)

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
58 of ANSYS, Inc. and its subsidiaries and affiliates.
 Model quality measures

The CV estimate is a random number that depends on the division into folds. Repeating cross-validation
multiple times using different splits into folds provides a better approximation to complete -fold
cross-validation (leave-one-out). Leave-one-out measure is almost unbiased, but for typical real world
datasets it has high variance, leading to unreliable estimates. Small datasets are simply not suitable
for CV estimates, since data distribution can change considerably after we separate out even a small
portion of data. In addition, the CV approach is usually too expensive. The question is whether the
advantages of CV (if any) are big enough to justify the computational cost of training multiple networks
rather than one.

3.2.4. Coefficient of Determination (R2)

Available in: DesignXplorer, LS-OPT, ModelCenter, optiSLang

The Coefficient of Determination (CoD) can be used to assess the approximation quality of a regression
model. This measure is defined as the relative amount of variation explained by the approximation
(Montgomery and Runger 2003 (p. 61)):

(3.70)

where is equivalent to the total variation of the output , represents the variation due to
the regression, and quantifies the unexplained variation,

(3.71)

If the CoD is close to one, the approximation represents the learning point values with small errors.
However, a polynomial model would fit exactly through the learning points, if their number is equi-
valent to the number of coefficients . In this case, the CoD would be equal to one, independent of
the true approximation quality. In order to penalize this over-fitting, based on the same approximation
model, the adjusted Coefficient of Determination was introduced (Montgomery and Runger 2003):
(3.72)

However, the over-estimation of the approximation quality can not be avoided completely.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 59
Metamodeling techniques

Figure 3.21: Subspace plot of the investigated nonlinear function (Equation 3.73 (p. 61)) and
convergence of the CoD measures with increasing number of learning points

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
60 of ANSYS, Inc. and its subsidiaries and affiliates.
 Model quality measures

In order to demonstrate this statement, an investigation of a nonlinear analytical function is performed.

The function of five independent and uniformly distributed input variables reads as follows:
(3.73)

where the contributions of the five inputs to the total variance are : 18.0%, : 30.6%, : 64.3%,
: 0.7%, : 0.2%. This means, that the three variables, , and , are the most important.

In Figure 3.21 (p. 60), the convergence of the standard CoD of linear and quadratic response surfaces
is shown, where a strong over-estimation of the approximation quality can be noticed, when the
number of samples is relatively small. Even the adjusted CoD shows a similar behavior. This fact limits
the CoD to cases where a large number of learning points compared to the number of polynomial
coefficients is available. However, in industrial applications this is often not the case. For other local
approximation models, like interpolating Kriging, this measure may be equal or close to one, however
the approximation quality may be poor.

References
Montgomery, D. C., and G. C. Runger. 2003. Applied Statistics and Probability for Engineers. Third. John
Wiley & Sons.

3.2.5. Coefficient of Prognosis

Available in: optiSLang

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 61
Metamodeling techniques

In (Most and Will 2008) a model independent measure to assess the model quality was proposed.
This measure is the Coefficient of Prognosis (CoP), which is defined as follows:

(3.74)

where is the sum of squared prediction errors. These errors are estimated based on cross
validation. In the cross validation procedure, the set of learning points is mapped to subsets. Then
the approximation model is built by removing subset from the learning points and approximating
the subset model output using the remaining point set. Based on definition (see Equation
3.68 (p. 58)), can be expressed as:
(3.75)
This means that the model quality is estimated only at those points which are not used to build the
approximation model. Since the prediction error is used instead of the fit, this approach applies to
regression and even interpolation models.

The evaluation of the cross validation subsets, which are usually between 5 and 10 sets, causes addi-
tional numerical effort in order to calculate the CoP. Nevertheless, for polynomial regression and
Moving Least Squares, this additional effort is still quite small since no complex training algorithm is
required. For other meta-modeling approaches such as neural networks, Kriging and even Support
Vector Regression, the time consuming training algorithm has to be performed for every subset
combination.

In Figure 3.22 (p. 62), the convergence of the CoP values for an MLS approximation of the nonlinear
coupled function given in Equation 3.73 (p. 61) is shown in comparison to the polynomial CoD. The
figure indicates that the CoP values are not over-estimating the approximation quality like the CoD
does for a small number of samples.

Figure 3.22: Convergence of the CoP measure by using MLS approximation compared to the
polynomial CoD measure

The influence radius of the MLS approximation is found by maximizing the CoP measure. As can be
seen in Figure 3.22 (p. 62), the convergence of the approximation quality is much better if only the
three most important variables are used in the approximation model.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
62 of ANSYS, Inc. and its subsidiaries and affiliates.
 Model Sensitivity measures

References
Most, T., and J. Will. 2008. "Metamodel of Optimal Prognosis - an Automatic Approach for Variable
Reduction and Optimal Metamodel Selection." In Proc. Weimarer Optimierungs- Und Stochastiktage 5.0,
Weimar, Germany, November 20-21, 2008.

3.3. Model Sensitivity measures

3.3.1. Variance based sensitivity analysis

Assuming a model with a scalar output as a function of a given set of random input parameters

(3.76)

the first order sensitivity index is defined as

(3.77)

where is the unconditional variance of the model output and is the variance of
caused by a variation of only.

Since first order sensitivity indices measure only the decoupled influence of each variable an extension
for higher order coupling terms is necessary. Therefore total effect sensitivity indices have been intro-
duced

(3.78)

where is the variance of caused by all model inputs without .

In order to estimate the first order and total sensitivity indices, a matrix combination approach is very
common (Saltelli et al. 2008). This approach calculates the conditional variance for each variable with
a new sampling set. In order to obtain a certain accuracy, this procedure requires often more than
1000 samples for each estimated conditional variance. Thus, for models with a large number of variables
and time consuming solver calls, this approach can not be applied efficiently.

References
Saltelli, A. et al. 2008. Global Sensitivity Analysis. The Primer. Chichester, England: John Wiley & Sons,
Ltd.

[Link]. ANOVA
Available in: LS-OPT

Since the number of regression coefficients determines the number of simulation runs, it is important
to remove those coefficients or variables which have small contributions to the design model. This
can be done by doing a preliminary study involving a design of experiments and regression analysis.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 63
Metamodeling techniques

The statistical results are used in an analysis of variance (ANOVA) to rank the variables for screening
purposes. The procedure requires a single iteration using polynomial regression, but results are
produced after every iteration of a normal optimization procedure. ANOVA is a regression based
sensitivity measure with

(3.79)

where is the linear approximation, the size of the design space of variable , and the
number of variables, Figure 3.23 (p. 64):

Figure 3.23: Definition of ANOVA value

[Link].1. The confidence interval of the regression coefficients

The confidence interval for the least squares estimators for is determined
by the inequality
(3.80)

where
(3.81)

is the regression coefficient and is an unbiased estimator of the variance given by

(3.82)

is the number of learning points and is the number of terms in the polynomial while
is the diagonal element of corresponding to and is Student"s -Distribution.
therefore represents the level of confidence that will be in the computed interval.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
64 of ANSYS, Inc. and its subsidiaries and affiliates.
 Model Sensitivity measures

Figure 3.24: ANOVA value with confidence interval

[Link].2. The significance of a regression coefficient

The contribution of a single regressor variable to the model can also be investigated. This is done
by means of the -test where is calculated to be

(3.83)

where =1 and the reduced model is the one in which the regressor variable in question has
been removed. Each of the terms represents the sum of squared residuals for the reduced and
complete models respectively.

It turns out that the computation can be done without analyzing a reduced model by computing

(3.84)

can be compared with the -statistic so that if , is non-zero with

confidence. The confidence level that is not zero can also be determined by computing the
for . The importance of is therefore estimated by both the magnitude of as
well as the level of confidence in a non-zero .

The significance of regressor variables may be represented by a bar chart of the magnitudes of
the coefficients with an error bar of length for each coefficient representing the
confidence interval for a given level of confidence . The relative bar lengths allow the analyst
to estimate the importance of the variables and terms to be included in the model while the error
bars represent the contribution to noise or poorness of fit by the variable.

All terms have been normalized to the size of the design space so that the choice of units becomes
irrelevant and a reasonable comparison can be made for variables of different kinds, e.g. sizing
and shape variables or different material constants.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 65
Metamodeling techniques

[Link]. Principal Component Analysis

Available in: LS-OPT, optiSLang

Principal component analysis (PCA) is a multivariate statistical technique using an orthogonal

transformation to convert a set of possibly correlated observations into a set of uncorrelated values
called principal components. The principal components are sorted in decreasing order according
to their variance, so that the first principal component accounts for the largest variability in the
data set (Abdi and Williams 2010).

Then, to analyze the effect of parameters on the data set, a method similar to ANOVA is applied
to the PCA results (Lamboni, Makowski, and Monod 2008).

The observations are represented by a matrix where is the number of experiments and
is the number of observations (number of discrete time values for histories, or number of points
for multi-responses).

The observation matrix is obtained using an approximation. For histories, the existing approximation
is used whereas for multi-responses a RBF approximation is used to get a full factorial representation
of the result.

The principal components are computed using the eigenvalue decomposition of the covariance
matrix (Missoum 2008),(Turk and Pentland 1991)
(3.85)

Where is the observation matrix centered and normalized, is the correlation matrix, are
the eigenvectors, and is the diagonal matrix of eigenvalues (in decreasing order) of the correl-
ation matrix. The Principal Components are computed as follow:
(3.86)

And is the inertia associated with the Principal Component. The sensitivity analysis of
the parameters is then performed by computing the variable sensitivity indices associated to the
principal components:

(3.87)

Where is the orthogonal projection matrix of (the variable values vector), is the sum
of squares associated to the variable for the principal component, and is the sensitivity
index of the variable to the principal component.

References
Abdi, H., and L. J. Williams. 2010. "Principal Component Analysis." Wiley Interdisciplinary Reviews:
Computational Statistics 2: 433–59.

Lamboni, M., D. Makowski, and D. Monod. 2008. "Multivariate Global Sensitivity Analysis for Discrete-
Time Models." {PhD} dissertationn, auto-saisine.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
66 of ANSYS, Inc. and its subsidiaries and affiliates.
 Model Sensitivity measures

Missoum, S. 2008. "Probabilistic Optimal Design in the Presence of Random Fields." Structural and
Multidisciplinary Optimization 35: 523–30.

Turk, M., and A. Pentland. 1991. "Eigenfaces for Recognition." Journal of Cognitive Neuroscience 3:
71–86.

[Link]. Coefficient of Correlation

Available in: optiSLang

The coefficient of correlation is the standardized covariance between two random variables and

(3.88)

where is the covariance and is the standard deviation. This quantity, known as the
linear correlation coefficient, measures the strength and the direction of a linear relationship between
two variables. It can be estimated from a given sampling set as follows:

(3.89)

where is the number of samples, and are the sample values, and and are the estimates
of the mean value and the standard deviation, respectively. The estimated correlation coefficient
becomes more inaccurate, as its value is closer to zero, which may cause a wrong deselection of
apparently unimportant variables.

If both variables have a strong positive correlation, the correlation coefficient is close to one. For
a strong negative correlation is close to minus one. The squared correlation coefficient can be
interpreted as the first order sensitivity index by assuming a linear dependence. The drawback of
the linear correlation coefficient is its assumption of just a linear dependence. Based on the estimated
coefficients only, it is not possible to decide on the validity of this assumption. Correlation coeffi-
cients, which assume a higher order dependence or use rank transformations solve this problem
only partially. Additionally, often interactions between the input variables are important. These in-
teractions can not be quantified with the linear and higher order correlation coefficients.

We can summarize that although the correlation coefficient can be simply estimated from a single
sampling set, it can only quantify first order effects with an assumed dependence without any
quality control of this assumption.

[Link]. Coefficient of Importance

Available in: optiSLang

The Coefficient of Importance (CoI) was developed to quantify the input variable importance by
using the CoD measure (Coefficient of Determination (R2) (p. 59)). Based on a polynomial model,
including all investigated variables, the CoI of a single variable with respect to the response
is defined as follows
(3.90)

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 67
Metamodeling techniques

where is the CoD of the full model including all terms of the variables in and is the
CoD of the reduced model, where all linear, quadratic and interactions terms belonging to are
removed from the polynomial basis. For both cases the same set of sampling points is used. If a
variable has low importance, its CoI is close to zero, since the full and the reduced polynomial re-
gression model have a similar quality. The CoI is equivalent to the explained variation with respect
to a single input variable, since the CoD quantifies the explained variation of the polynomial ap-
proximation. Thus it is an estimate of the total effect sensitivity measure given in Equation
3.78 (p. 63). If the polynomial model contains important interaction terms, the sum of the CoI values
should be larger than the CoD of the full model.

Since it is based on the CoD, the CoI is also limited to polynomial models. If the total explained
variation is over-estimated by the CoD, the CoI may also give a wrong estimate of the variance
contribution of the single variables. However, in contrast to the Coefficient of Correlation, the CoI
can handle linear and quadratic dependencies including input variable interactions. Furthermore,
an assessment of the suitability of the polynomial basis is possible. Nevertheless, an estimate of
the CoI values using a full quadratic polynomial is often not possible because of the required large
number of samples for high dimensional problems.

3.4. MOP principle and competition

Available in: optiSLang

As demonstrated in section Section 3.2.5 (p. 61), the prediction quality of an approximation model may
be improved if unimportant variables are removed from the model. This idea is adopted in the
Metamodel of Optimal Prognosis (MOP) proposed in (Most and Will 2008) which is based on the search
for the optimal input variable set and the most appropriate approximation model. Due to the model
independence and objectivity of the CoP measure, it is well suited to compare the different models in
the different subspaces.

Figure 3.25: CoP values of different input variable combinations and approximation methods
obtained with the analytical nonlinear function

In Figure 3.25 (p. 68), the CoP values of all possible subspaces and all possible approximation models
are shown for the analytical nonlinear function example (Equation 3.73 (p. 61)). The figure indicates
that there exists an optimal compromise between the available information, the learning points and

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
68 of ANSYS, Inc. and its subsidiaries and affiliates.
 MOP principle and competition

the model complexity i.e. the number of input variables. The MLS approximation by using only the
three major important variables has a significantly higher CoP value than other combinations. However
for more complex applications with many input variables it is necessary to test a huge number of ap-
proximation models. In order to decrease this effort, in (Most and Will 2008) advanced filter technologies
are proposed, which reduce the number of necessary model tests. Nevertheless, a large number of inputs
requires a very fast and reliable construction of the approximation model. For this reason polynomials
and MLS are preferred due to their fast evaluation.

As a result of the MOP, we obtain an approximation model which contains the most important variables.
Based on this meta-model, the total effect sensitivity indices, proposed in section Section 3.3.1 (p. 63),
are used to quantify the variable importance. The variance contribution of a single input variable is
quantified by the product of the CoP and the total effect sensitivity index estimated from the approx-
imation model
(3.91)

Since interactions between the input variables can be represented by the MOP approach, they are
considered automatically in the sensitivity indices. If the sum of the single indices is significantly larger
as the total CoP value, such interaction terms have significant importance.

Additionally to the quantification of the variable importance, the MOP can be used to visualize the de-
pendencies in 2D and 3D subspaces. This helps the designer to understand and to verify the solver
model. In Figure 3.26 (p. 69) two subspace plots are shown for the MOP of the analytical test function
(Equation 3.73 (p. 61)). In the - and - subspace plots, the sinusoidal function behavior and
the coupling term can be observed.

Figure 3.26: - and - subspace plots of the MOP of the nonlinear analytical function given
in Equation 3.73 (p. 61)

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 69
Metamodeling techniques

Additional parametric studies, such as global optimization can also be directly performed on the MOP.
Nevertheless, a single solver run should be used to verify the final result of such a parametric study or
optimization.

References
Most, T., and J. Will. 2008. "Metamodel of Optimal Prognosis - an Automatic Approach for Variable Re-
duction and Optimal Metamodel Selection." In Proc. Weimarer Optimierungs- Und Stochastiktage 5.0,
Weimar, Germany, November 20-21, 2008.

3.4.1. MOP residual plot

Available in: optiSLang

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
70 of ANSYS, Inc. and its subsidiaries and affiliates.
 MOP principle and competition

Figure 3.27: Residual plot of the MOP-post-processing: original vs. the approximated values
(left) and the sample CoP indicated for each design (right)

In the MOP post-processing the residual plot is a simple way to investigate the prediction accuracy
of the individual designs. In this plot, both the approximated values ("Fitting" layer) and the test values
from cross validation ("Prediction" layer) are plotted against original data values. For a perfect model,
the points will collapse onto a straight identity line through origin (black line). Additionally, the
global maximum, mean and root mean squared error measures as well as the Coefficient of Determ-
ination and Coefficient of Prognosis are indicated for the fitting and the cross-validation residuals.
An additional error measure is provided in the residual plot, the Sample Coefficient of Prognosis,
which indicates the contribution of each sample to the global CoP. The calculation for a specific
sample reads as follows:

(3.92)

The mean of the Sample CoPs of all samples is equal to the global CoP, if each individual Sample CoP
is larger as zero.

3.4.2. Local MOP error measures

Available in: optiSLang

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 71
Metamodeling techniques

Figure 3.28: MOP-post-processing: local CoP calculated as weighted averaging of the sample
CoPs

Figure 3.29: Estimated support density of a design set with almost uniform distribution (left)
and a pure random distribution with regions with high and low density (right)

In the 2D and 3D response surface plots of the MOP post-processing local approximation errors are
available as indicated in Figure 3.28 (p. 72). These local error measures are directly obtained from
the cross validation residuals of the individual samples by using a local averaging scheme similar to
the Moving Least Squares approximation. The locally weighted root mean squared error and the
local CoP read as follows

(3.93)

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
72 of ANSYS, Inc. and its subsidiaries and affiliates.
 Adaptive metamodeling

(3.94)

An additional important measure in the meta-model post-processing is the support density. This
measure gives an idea how well the support points used for the meta-model generation are distributed
in the design space. As shown in Figure 3.29 (p. 72) for an almost uniformly distributed sample set,
the support density is close to one. In a pure random set, this density measures indicates clusters
with higher values as two and indicates a very low density with values below 0.1.

3.5. Adaptive metamodeling

3.5.1. Adaptive Metamodel of Optimal Prognosis (AMOP)

Available in: optiSLang

Figure 3.30: Adaptive MOP: approximation history plot with the convergence of the global and
minimum local CoP values for each response

In the Adaptive MOP approach, the initial meta-model trained on a given data set, is improved step-
by step with new samples in those regions, where the benefit is optimal for the later application.

In the global refinement, the initial DoE is improved by an additional space-filling or correlation op-
timized Latin Hypercube Sampling in a global update scheme (see section Latin Hypercube
Sampling (p. 12)). Since the two advanced LHS DoE-schemes consider the existing designs, either
the spurious correlations or the space filling property are optimized. If the global CoPs for all selected
responses reach the specified target CoP, the iteration is stopped. The user can specify the number
of samples for the initial iteration or can consider existing designs from previous analysis systems.
Additionally, the user has to define the number of samples and the maximum number of the following
iterations.

As local refinement strategies, a pure density based method, a quality based approach and a criteria
refinement are available and can be combined with each other by user-defined weighting factors. In
the pure density and quality-based strategies, new samples are generated in regions with low density

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 73
Metamodeling techniques

or low local approximation quality. The local approximation error is directly estimated based on the
cross-validation residuals of the initial data points as shown in Figure 3.28 (p. 72). If the minimum
sample CoP of all selected responses reaches the target CoP, the iteration is stopped. In Figure
3.30 (p. 73) an example of the convergence of the global and local CoP values is shown. In Figure
3.31 (p. 74) the corresponding update scheme is illustrated in the design space of the simple 2D ex-
ample. It can be seen in the figure, that the new samples are mainly placed in the region with the
lowest local CoPs which improves the local approximation quality significantly.

Figure 3.31: Adaptive MOP: Local refinement considering the local approximation quality

Figure 3.32: Adaptive MOP: Local refinement considering constraints of the response values

Another update strategy uses the constraint definitions for the investigated response values with
specified value ranges for the model improvement, as shown in Figure 3.32 (p. 74). In this refinement
strategy, the samples are placed in regions, where the constraints are assumed to be fulfilled using
a space-filling in-fill criterion.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
74 of ANSYS, Inc. and its subsidiaries and affiliates.
 Adaptive metamodeling

The local refinement strategy of the AMOP can also be used for single and multi-objective optimization
and is introduced in the optimization section.

3.5.2. Genetic Aggregation Response Surface Refinement

Available in: DesignXplorer Gaussian Process regression methods, such as Kriging, provide a prediction
distribution. The genetic aggregation can contain several surrogate models that do not afford a local
prediction distribution. Ben Salem et al. (Ben Salem et al. 2017) introduced a universal prediction (UP)
for all surrogate models. This distribution relies on cross-validation sub-model predictions. Based on
this empirical distribution, they define a variance (UP variance) and an adaptive sampling technique
for global refinement (UP-SMART).

Figure 3.33: Illustration of the UP distribution for a kriging surrogate (left). Dashed lines: CV
sub-models predictions, solid red line: master model prediction, horizontal bars: local UP
distribution at = -1.8 and = 0.2, black squares: design points.

Figure 3.34: Uncertainty quantification based on the UP distribution for a kriging surrogate.
Blue solid line: master model prediction , light blue area: region delimited by

The Genetic Aggregation refinement is a hybrid variant of UP-SMART consisting in adding at step
a point:
(3.95)

Where:

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 75
Metamodeling techniques

• is the local UP variance at point et the step

• and denotes a given distance on

• is the sample set used to build the metamodel

• is the parameter space

• is the distance penalization

Remark 1. The predicted error exposed in DesignXplorer corresponds to .

The first refinement point selected is the worst point, having the maximum absolute predicted error
(APE) in regard to the expected tolerance.

Figure 3.35 (p. 76) represents the absolute predicted error of a response surface with only one input
parameter.

Figure 3.35: Absolute predicted error in a 1-dimensional problem. Circles: design points used
to build the metamodel, diamond: 1st refinement point

For refinement points to submit simultaneously, the generation of the refinement point depends
on the ) first pending refinement points, with . The response surface is updated only when
all refinement points have been generated.

Figure 3.36: Absolute weighted predicted error in a 1-dimensional problem. Circles: design
points used to build the metamodel, diamond: 1st refinement point, triangle: 2nd refinement
point

While the previous figure showns that the first refinement point is based on the APE, in the next figure,
refinement points are based on the absolute weighted predicted error (AWPE), which takes into account
the influence of the pending refinement points on the response surface. AWPE is a transformation
of APE that reduces the predicted error around the pending refinement point to favor the domain

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
76 of ANSYS, Inc. and its subsidiaries and affiliates.
 Metamodels Glossary

with a high predicted error and a low density of pending refinement points. The refinement point,
with , is selected as the worst point with the maximum AWPE.

References
Ben Salem, Malek, Olivier Roustant, Fabrice Gamboa, and Lionel Tomaso. 2017. "Universal Prediction
Distribution for Surrogate Models." SIAM/ASA Journal on Uncertainty Quantification 5 (1): 1086–1109.

3.6. Metamodels Glossary

ANOVA. Analysis of variance. Used to perform variable screening by identifying insignificant variables.
Variable regression coefficients are ranked based on their significance as obtained through a partial F-
test. (See also variable screening).

Bias error. The total error - the difference between the exact and computed response - is composed
of a random and a bias component. The bias component is a systematic deviation between the chosen
model (approximation type) and the exact response of the structure (FEA analysis is usually considered
to be the exact response). Also known as the modeling error. (See also random error).

Confidence interval. The interval in which a parameter may occur with a specified level of confidence.
Computed using Student"s t-test. Typically applied to accompany the significance of a variable in the
form of an error bar.

Dependent. A function which is dependent on variables. Dependent variable.

Design of Experiments. See experimental design.

Design parameter. See design variable.

Design space. A region in the -dimensional space of the design variables ( through ) to which
the design is limited. The design space is specified by upper and lower bounds on the design variables.
Response variables can also be used to bound the design space.

Design surface. The response variable as a function of the design variables, used to construct the for-
mulation of a design problem. (See also response surface).

Design variable. An independent design parameter which is allowed to vary in order to change the
design. Symbolized by ( or (vector containing several design variables)).

DFFN. Deep Feed Forward Network. A Metamodeling technique

DIM-GP. Deep Infinite Mixture Gaussian Process. A Metamodeling technique.

DOE. Design of Experiments. See experimental design.

Experimental Design. The selection of designs to enable the construction of a design response surface.
Sometimes referred to as the Point Selection Scheme.

Feedforward Neural Network. See Neural Network.

Function. A mathematical expression for a response variable in terms of design variables. Often used
interchangeably with "response". Symbolized by f.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 77
Metamodeling techniques

Function evaluation. Using a solver to analyze a single design and produce a result. See Simulation.

FMU. Functional Mock-Up Unit. A simulation model or program which implements the FMI standard.

FMI. Functional Mock-up Interface. The FMI standard is a tool independent standard to support both
model exchange and co-simulation of dynamic models using a combination of XML-files, C-header files,
C-code or binaries.

Global approximation. A design function which is representative of the entire design space.

GARS. Genetic Aggregation Response Surface. A Metamodeling technique.

Iteration. A cycle involving an experimental design, function evaluations of the designs, approximation
and optimization of the approximate problem.

Kriging. A Metamodeling technique using Bayesian regression.

Metamodeling. The construction of surrogate design models such as polynomial response surfaces,
Artificial Neural Networks or Kriging surfaces from simulations at a set of design points.

Modeling error. See bias error.

MLS. Moving Least Squares. A Metamodeling technique.

MOP. Metamodel of Optimal Prognosis. General framework for an automatic model testing and selection
in optiSLang.

MSE. Mean Squared Error.

Neural network approximation. The use of trained feedforward neural networks to perform non-linear
regression, thereby constructing a non-linear metamodels (see metamodeling).

Noise. See random error.

Radial basis function network. The use of radial basis functions (RBFs) to approximate response
functions. This is a global approximation method.

Random error. The total error - the difference between the exact and computed response - is composed
of a random and a bias component. The random component is, as the name implies, a random deviation
from the nominal value of the exact response, often assumed to be normally distributed around the
nominal value. (See also bias error).

Region of interest. A sub-region of the design space. Usually defined by a mid-point design and a
range of each design variable. Usually dynamic.

Residual. The difference between the computed response (using simulation) and the predicted response
(using a response surface).

Response. A numerical indicator of the performance of the design. A function of the design variables
approximated using a metamodel which can be used for optimization. Symbolized by f. Collected over
all design iterations for plotting.

Response quantity. See response.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
78 of ANSYS, Inc. and its subsidiaries and affiliates.
 Metamodels Glossary

Response Surface. A mathematical expression which relates the response variables to the design
parameters. Typically computed using statistical methods.

Response variable. A response which is substituted as an input variable in a downstream stage of a

multi-stage process. Response-expressions are allowed.

Result. A numerical indicator of the performance of the design. A result is not associated with a
metamodel, but is typically used for intermediate calculations in metamodel-based analysis.

RBF. Radial Basis Function. RBF"s are used as basis functions for metamodels (see also metamodeling).
These functions are typically Gaussian.

RSM. Response Surface Methodology.

Saturated design. An experimental design in which the number of points equals the number of unknown
coefficients of the approximation. For a saturated design no test can be made for the lack of fit.

Sensitivity. See Design sensitivity.

Simulation. The analysis of a physical process or entity in order to compute useful responses. See
Function evaluation.

Slack variable. The variable which is minimized to find a feasible solution to an optimization problem,
e.g. in: min subject to . See Strictness.

Solver. A computational tool used to analyze a structure or fluid using a mathematical model. The ex-
ecutable software for such a tool.

Stochastic. Involving or containing random variables. Involving probability or chance.

Strictness. A number between 0 and 1 which signifies the strictness with which a design constraint
must be treated. A zero value implies that the constraint may be violated. If a feasible design is possible
all constraints will be satisfied. Used in the design formulation to minimize constraint violations. See
Slack variable.

Subregion. See region of interest.

SVR. Support Vector Regression. A Metamodeling technique.

Variable screening. Method to remove insignificant variables from the design optimization process
based on a ranking of regression coefficients using analysis of variance (ANOVA). (See also ANOVA).

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 79
 Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
80 of ANSYS, Inc. and its subsidiaries and affiliates.
Chapter 4: Optimization methods
The following topics are available:
4.1. Optimization setup
4.2. Single-objective optimization
4.3. Multi-objective optimization
4.4. Optimization Glossary

4.1. Optimization setup

Optimization can be defined as a procedure for achieving the best outcome of a given operation while
satisfying certain restrictions. This objective has always been central to the design process, but is now
assuming greater significance than ever because of the maturity of mathematical and computational
tools available for design.

In single and multi-objective optimization the optimization task has to be defined by means of the
objective functions:
(4.1)

which can be implicit functions of the input variables and of scalar and non-scalar response values
. In optiSLang the user can define minimization and maximization task with scalar objective functions,
where the full calculator functionality can be applied. For non-scalar response values such as signal
outputs, scalarization functions are available to derive scalar measures for the definition of the objectives
such as integral values, deviation between two signals and many more.

Figure 4.1: Definition of optimization parameter properties in optiSLang

In unconstrained optimization problems only the bounds or values of the design variables limit the
optimization space. The optimizer searches between these limits for the minimum value of the objective
function . For maximization problems the sign of the user-defined objective function is inverted
in optiSLang to get a minimization task. The design variables can be defined as continuous variables
with a lower and upper bound or as discrete variables which assume several discrete values. In Figure
4.1 (p. 81) the definition of different optimization parameters types for optiSLang is shown. The discrete
parameters are distinguished as discrete by value, which use the real data axis, ordinal discrete, which

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 81
Optimization methods

use the ordered ordinal axis, and nominal discrete, which may contain completely unordered states
having real, integer or string value types.

In engineering problems often additional restrictions have to be fulfilled by the optimal design. With
help of equality and inequality constraints:

(4.2)

such restrictions can be formulated. The constraint functions can represent limitations depending only
on the input variables but also depending on all available model responses and any mathematical
combination of both.

Figure 4.2: Definition of objective functions and optimization constraints in optiSLang

Equality constraints are very difficult to treat. Therefore, it is not allowed to define equality constraints
in optiSLang. However, equality constraints fulfilled with a given tolerance can be formulated easily as
inequality constraints and in that way they can be handled by all optimization methods available in
optiSLang. In order to get an optimal convergence of the optimizer and sufficiently accurate fulfillment
of the constraint conditions, it is very useful to scale the objective and the constraint functions to be
in a similar range. This can be realized e. g. by using the results of the design exploration.

In Figure 4.2 (p. 82) the typical setup of the optimization criteria, objectives functions and constraints,
is shown. The objective functions can be defined as maximization or minimization goals, whereas the
constraint function contain a left and a right functional term within a greater equal or less equal
definition.

Table 4.1 (p. 82) describes the properties and a suggested application for each optimizer introduced in
this documentation.

Table 4.1: Available optimization algorithms with application for single and multi-objective
optimization and possible discrete inputs variables

Optimizer Objectives Discrete

inputs

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
82 of ANSYS, Inc. and its subsidiaries and affiliates.
 Optimization setup

Single Multi by value ordinal nominal

Hybrid methods
One Click Optimization yes yes yes yes no
Gradient-based methods
Non-Linear Programming (NLPQLP) yes no no no no
Mixed Integer SQP yes no yes yes no
Leapfrog Optimization (LFOPC) yes no (yes) (yes) no
DAKOTA Coliny Solis-Wets yes no no no no
DAKOTA Coliny CONMIN yes no no no no
DAKOTA Coliny FDN yes no no no no
DAKOTA Coliny Quasi Newton (QN) yes no no no no
DAKOTA Coliny PR yes no no no no
Pattern search methods
Downhill simplex yes no (yes) (yes) no
Hooke-Jeevs Pattern Search yes no no no no
DAKOTA OPT++ Parallel Direct Search yes no no no no
DAKOTA APPS yes no no no no
DAKOTA Coliny DIRECT yes no no no no
DAKOTA Coliny Pattern Search yes no no no no
Response surface based methods
Adaptive Response Surface Method (ARSM) yes no yes yes no
Adaptive MOP (AMOP) yes yes yes yes no
Adaptive Single-Objective (ASO) yes no yes yes no
Adaptive Multi-Objective (AMO) no yes yes yes no
Efficient Global Optimization (UP-EGO) yes no yes yes no
Bayesian Optimization (PI-BO) yes no yes yes (yes)
Nature-inspired methods
Evolutionary Algorithm (EA) yes yes yes yes yes
Darwin algorithm yes yes yes yes yes
DAKOTA Coliny Evolutionary Algorithm yes no yes yes yes
EVOLVE yes no yes yes yes
Covariance Matrix Adaptation (CMA) yes no yes yes (yes)
Particle Swarm Optimization (PSO) yes yes yes yes (yes)
Stochastic Design Improvement (SDI) yes no yes yes yes
Adaptive Simulated Annealing yes no yes yes (yes)
Differential Evolution yes no yes yes (yes)

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 83
Optimization methods

4.2. Single-objective optimization

In single-objective optimization problems the optimization task is formulated just by a single scalar-
valued objective function:
(4.3)

Figure 4.3: Recommended flowchart for single-objective optimization

In Figure 4.3 (p. 84) the recommended flow of single-objective optimization procedure is shown: after
the definition of the design variables and objective and constraint functions the design space is explored
by sensitivity analysis. The obtained variable sensitivities may help to reduce the number of design
variables. The best designs found in the sensitivity analysis could be used as start designs for the fol-
lowing optimization procedure which finally will determine an optimal design.

4.2.1. Accompanying example: optimization of a damped oscillator

The principles of the optimization methods described in the following sections are demonstrated by
a simple accompanying example. For this purpose the maximum amplitude of a single degree-of-
freedom damped oscillator shown in Figure 4.4 (p. 84) shall be minimized.

Figure 4.4: Damped oscillator: system properties and oscillation behavior

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
84 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

The equation of motion can be formulated depending on the mass , the spring stiffness and
damping ratio as follows:
(4.4)

where is the undamped eigen-frequency. Assuming an initial kinetic energy and an

initial position , the time-dependent displacement function can be derived as follows:
(4.5)

where is the damped eigen-frequency. The optimization goal is to minimize the max-
imum amplitude after 5 seconds by assuming and as design parameters and and as
constants

(4.6)

As optimization constraint, the damped eigen-frequency shall be limited: . If the maximum

amplitude after 5 seconds is estimated from the envelope curve shown in Figure 4.4 (p. 84), the ob-
jective function is a smooth and well defined function of the design parameters as shown in Figure
4.5 (p. 86).

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 85
Optimization methods

Figure 4.5: Objective function of the damped oscillator obtained by using an estimate of the
maximum amplitude from the envelope curve (left) and by using a coarse time discretization
of the displacement curve (right)

The constraint condition is indicated in the figure by the function . If the maximum amplitude
is obtained by a time discretization of Equation 4.5 (p. 85) using a coarse time interval, the resulting
objective function may contain local oscillations, which may be interpreted as solver noise. In order
to demonstrate the behavior of the different optimizers in the presence of small solver noise, a time
step of 0.1 seconds is chosen, which leads to a slightly noisy objective function shown additionally
in Figure 4.5 (p. 86).

4.2.2. Gradient-based methods

Available in: DesignXplorer, LS-OPT, ModelCenter, optiSLang

Gradient-based optimization methods use local derivatives of the objective function to find the next
local optimum. The minimum of a convex function can be determined by searching for the point,
where the first derivative is equal to zero as shown in Figure 4.6 (p. 87). If the second derivatives is
known, the objective function can be approximated by a second order Taylor series. This procedure
is called Nonlinear Programming (NLP). Optimization methods using first and second order derivatives
are known as Newton methods (Kelley 1999).

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
86 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

Figure 4.6: Iterative search of a gradient-based method using a quadratic approximation of the
objective function

If the CAE solver is treated as black box, the derivatives have to be calculated by numerical estimates.
Since second order derivatives would require a large number of solver calls, full Newton methods are
not applicable for complex optimization problems. For this reason quasi Newton methods have been
developed, where the second order derivatives are approximated from the first order derivatives of
previous iteration steps. Gradient-based methods distinguish each other mainly in the way how the
second order derivatives are estimated and how additional constraint equations are considered. A
good overview over different methods in given in (Kelley 1999).

References
Kelley, C. T. 1999. Iterative Methods for Optimization. Philadelphia: Siam, Society for Industrial; Applied
Mathematics.

[Link]. Non-Linear Programming by Quadratic Lagrangian (NLPQLP)

Available in: DesignXplorer, ModelCenter, optiSLang

Keywords— Mathematical programming, SQP

NLPQLP solves smooth nonlinear programming problems by a sequential quadratic programming

(SQP) algorithm (Schittkowski 1986). SQP is a class of algorithms for solving non-linear optimization
problems (NLP). The theoretical background is described e. g. in Stoer (Stoer 1985), Fletcher

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 87
Optimization methods

(Fletcher 1981), or Gill et al.(GILL, MURRAY, and WRIGHT 1981). First-order derivatives of the objective
and constraint functions are estimated from central differences or by one-sided numerical derivatives
using a specified differentiation interval.

The NLPQL approach features efficient constraint handling using quadratic Lagrangian. Starting
from a given start point, the method searches for the next local optimum and converges if the es-
timated gradients are below a specified tolerance. Since it is a local optimization method, it is re-
commended to use the best design of a global sensitivity analysis as start point in order to find
the global optimum. The NLPQL is very efficient in low dimensions with up to 20 design variables.
For higher dimensional problems the computation of the numerical derivatives becomes more and
more expensive and other methods are more efficient. In the presence of noisy model responses
the differentiation interval plays a crucial role. If it is taken too small, the estimated gradient is dis-
torted heavily by the solver noise and the NLPQL runs into a wrong direction. If the solver noise is
not dominating the functional trend an increased differentiation interval could lead to a good
convergence of the optimization procedure.

Figure 4.7: Convergence of the NLPQL optimizer for the oscillator optimization problem by
using the smooth objective function based on the envelope estimate (left) and by using the
noisy objective function from the time discretization (right)

Example: In Figure 4.7 (p. 88) the convergence of the NLPQL approach is shown for the optimization
of the damped oscillator: in the first investigation, the smooth objective function obtained from
the envelope curve estimate is considered. In this case the optimizer converges in three iteration
steps to the true optimum ( , ). If the noisy objective function is investigated
with a relatively small differentiation interval ( of the design space) the optimizer runs into the
wrong direction and will not converge to the true optimum. Nevertheless, by increasing the differ-
entiation interval the convergence to the true optimum can be achieved for this example. Limita-
tions: NLPQL does not support discrete input parameters.

References
Fletcher, R_. 1981. "Practical Methods of Optimization: Vol. 2: Constrained Optimization." JOHN
WILEY & SONS, INC., ONE WILEY DR., SOMERSET, N. J. 08873, 1981, 224.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
88 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

GILL, PE, W MURRAY, and MH WRIGHT. 1981. "Practical Optimization(book)." London and New York,
Academic Press, 1981. 415 p.

Schittkowski, K. 1986. "NLPQL: A Fortran Subroutine for Solving Constrained Nonlinear Programming
Problems." Annals of Operations Research 5: 485–500.

Stoer, J. 1985. "Foundations of Recursive Quadratic Programming Methods for Solving Nonlinear
Programs." Computational Mathematical Programming 15.

[Link]. Mixed-Integer Sequantial Quadratic Programming (MISQP)

Available in: DesignXplorer, optiSLang

Keywords— Mathematical programming, SQP

MISQP (Mixed-Integer Sequential Quadratic Programming) is a mathematical optimization algorithm

developed by Oliver Exler, Thomas Lehmann and Klaus Schittkowski (Oliver Exler and Schittkowski
2007),(O. Exler, Lehmann, and Schittkowski 2012). This method solves mixed-integer non-linear
programming problems (MINLP) by a modified sequential quadratic programming (SQP) method.
After linearizing constraints and constructing a quadratic approximation of the Lagrangian function,
mixed-integer quadratic programs are successively generated and solved by an efficient branch-
and-cut method. The algorithm is stabilized by a trust region method as originally proposed by
Yuan for continuous programs (Yuan 1995). Second order corrections are retained. The Hessian of
the Lagrangian function is approximated by BFGS updates subject to the continuous and integer
variables. MISQP is able to solve also non-convex nonlinear mixed-integer programs.

MISQP solves MINLP problems of the form:

(4.7)

The symbols and denote the vectors of the continuous and integer variables, respectively. It
is assumed that problem functions and are continuously differentiable
subject to all .

Limitations: It is not assumed that integer variables can be relaxed. In other words, problem
functions are evaluated only at integer points and never at any fractional values in between.

References
Exler, O., T. Lehmann, and K. Schittkowski. 2012. "MISQP: A Fortran Subroutine of a Trust Region
SQP Algorithm for Mixed-Integer Nonlinear Programming - Users Guide." Department of Computer
Science, University of Bayreuth.

Exler, Oliver, and Klaus Schittkowski. 2007. "A Trust Region SQP Algorithm for Mixed-Integer Nonlinear
Programming." Optimization Letter 1 (2): 269–80.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 89
Optimization methods

Yuan, Ya-xiang. 1995. "On the Convergence of a New Trust Region Algorithm." Numerische Mathematik
70 (4): 515–39.

[Link]. Leapfrog optimizer for constrained minimization (LFOPC)

Available in: LS-OPT

Keywords— Mathematical programming, Direct Search

The LFOPC algorithm solves constrained optimization problems (J. Snyman 2000). It is a gradient
method that generates a dynamic trajectory path, from any given starting point, towards a local
optimum. This method differs conceptually from other gradient methods, such as SQP, in that no
explicit line searches are performed. The original leap-frog method (J. A. Snyman 1983) for uncon-
strained minimization problems seeks the minimum of a function of n variables by considering the
associated dynamic problem of a particle of unit mass in an n-dimensional conservative force field,
in which the potential energy of the particle at point x(t) at time t is taken to be the function
to be minimized. The solution to the constrained problem may be approximated by applying the
unconstrained minimization algorithm to a penalty function formulation of the original algorithm.
The LFOPC algorithm uses a penalty function formulation to incorporate constraints into the optim-
ization problem. This implies that when constraints are violated (active), the violation is magnified
and added to an augmented objective function, which is solved by the gradient-based dynamic
leap-frog method (LFOP).

The algorithm uses three phases: Phase 0, Phase 1 and Phase 2. In Phase 0, the active constraints
are introduced as mild penalties through the pre-multiplication of a moderate penalty parameter
value. This allows for the solution of the penalty function formulation where the violation of the
(active) constraints are pre-multiplied by the penalty value and added to the objective function in
the minimization process. After the solution of Phase 0 through the leap-frog dynamic trajectory
method, some violations of the constraints are inevitable because of the moderate penalty. In the
subsequent Phase 1, the penalty parameter is increased to more strictly penalize violations of the
remaining active constraints. Finally, and only if the number of active constraints exceed the number
of design variables, a compromised solution is found to the optimization problem in Phase 2. Oth-
erwise, the solution terminates having reached convergence in Phase 1.

The values of the responses are scaled with the values at the initial design. The variables are scaled
internally by scaling the design space to [0; 1] interval. The default parameters in LFOPC should
therefore be adequate.

The optimization is terminated when either of the convergence criteria becomes active that is when
(4.8)

where refers to the vector of design variables and is the size of the design space. Or when

(4.9)

where denotes the value of the objective function, ( ) and ( -1) refer to two successive iteration
numbers.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
90 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

In the case of an infeasible optimization problem, the solver will find the most feasible design
within the given region of interest bounded by the simple upper and lower bounds. A global
solution is attempted by multiple start designs among the initial designs generated by the DoE
(see in Design of Experiments (p. 3)).

References
Snyman, JA. 2000. "The LFOPC Leap-Frog Algorithm for Constrained Optimization." Computers &
Mathematics with Applications 40 (8-9): 1085–96.

Snyman, Johannes Arnoldus. 1983. "An Improved Version of the Original Leap-Frog Dynamic
Method for Unconstrained Minimization: LFOP1 (b)." Applied Mathematical Modelling 7 (3): 216–18.

[Link]. DAKOTA Coliny Solis-Wets

Available in: ModelCenter

Keywords— Mathematical programming

The Solis-Wets is a heuristic local search algorithm for continuous design variables. Solis Wets
generates trial points using a multivariate normal distribution, and unsuccessful trial points are re-
flected about the current point to find a descent direction. This is a non-derivative optimization
algorithm and constraint violations are handled by a simple penalty function.

[Link]. DAKOTA CONMIN

Available in: ModelCenter

Keywords— Mathematical programming

CONMIN is a gradient based optimizer that uses the method of feasible direction to solve constrained
problems. CONMIN was a predecessor to the DOT optimizer.

[Link]. DAKOTA OPT++ Finite Differences Newton (FDN)

Available in: ModelCenter

Keywords— Mathematical programming

Newton method assumes that the objective function can be approximated as a quadratic function
around the local optimum. This algorithm uses a finite difference scheme to compute gradients
and the 2nd derivative matrix (Hessian). Constraints are handled by an interior-point penalty function.

[Link]. DAKOTA OPT++ Quasi Newton (QN)

Available in: ModelCenter

Keywords— Mathematical programming

Newton method assumes that the objective function can be approximated as a quadratic function
around the local optimum. This algorithm uses approximate Hessian (2nd derivative matrix) that is
updated as iteration progresses. Constraints are handled by an interior-point penalty function.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 91
Optimization methods

[Link]. DAKOTA OPT++ Polak-Ribiere Conjugate Gradient (PR)

Available in: ModelCenter

Keywords— Mathematical programming

DAKOTA OPT++ Conjugate Gradient is an unconstrained algorithm. The algorithm does not support
bounds on variables. Conjugate method assumes that the objective function can be approximated
as a quadratic function around the local optimum. This method uses the quadratic assumption to
calculate the search direction calculated

4.2.3. Pattern search methods

Available in: ModelCenter, optiSLang

Pattern search optimization is a type of optimization method that does not require gradient information
to locate the minimum. Instead, it uses a pattern of points to explore the search space and iteratively
refines the search by successively dividing the search space into smaller regions, or simplices, as it
approaches the minimum. This method is particularly useful when the optimization problem does
not have derivative information or when the problem is noisy or discontinuous.

[Link]. Downhill simplex

Available in: ModelCenter, optiSLang

Keywords— Mathematical programming, Direct Search

The downhill simplex algorithm, also known as the Nelder-Mead method (Nelder and Mead 1965),
is an iterative optimization algorithm used to find the minimum (or maximum) of a multivariate
function. It is a direct search method that does not require the derivative information of the function.
The algorithm starts with an initial simplex, which is a geometrical shape defined by a set of points
in the parameter space. In the case of a two-dimensional problem, the simplex is a triangle, while
in higher dimensions, it takes the form of a tetrahedron or higher-dimensional polytope. Each point
of the simplex represents a possible solution to the optimization problem. The algorithm iteratively
explores the parameter space by performing a series of transformations on the simplex. At each
iteration, it evaluates the objective function at the vertices of the simplex and uses this information
to determine which transformations to apply.

Step-by-step explanation of the downhill simplex algorithm:

1. Initialization: Choose an initial simplex, which is a set of points in the search space. The simplex
is typically an +1 dimensional shape, where is the number of variables in the function. Each
point in the simplex represents a solution to the optimization problem.

2. Order the points in the simplex according to their function values, with the highest value at
one vertex and the lowest value at another vertex.

3. Reflection (see in Figure 4.8 (p. 93)): Compute the reflection point ( ) by reflecting the worst
point ( , highest function value) across the centroid ( ) of the remaining points. The centroid
( ) is the average of the n best points.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
92 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

4. Evaluate the function at the reflection point. If the reflected point ( ) has a better function
value than the second worst point (but not better than the best point), proceed to step 6.

5. Expansion (see in Figure 4.8 (p. 93)): If the reflected point ( ) has a better function value than
the best point, compute an expansion point ( ) by extending the reflection point even further
in the same direction. Evaluate the function at the expansion point ( ). If the expansion point
( ) has a better function value than the reflected point ( ), replace the worst point ( ) in the
simplex with the expansion point ( ). Otherwise, replace the worst point ( ) with the reflected
point ( ).

6. Contraction (see in Figure 4.8 (p. 93)): If the reflected point ( ) does not have a better function
value than the second worst point, compute a contraction point ( ) by contracting the simplex
towards the best point. Evaluate the function at the contraction point ( ). If the contraction
point ( ) has a better function value than the worst point ( ), replace the worst point with
the contraction point. Otherwise, proceed to step 7.

7. Shrink (see in Figure 4.8 (p. 93)): If none of the above steps result in a better function value
than the worst point, perform a shrink operation. This involves contracting the entire simplex
towards the best point ( ). Each new point ( ) in the shrunk simplex is computed by taking
the average of the best point and one of the other points.

The steps 2-7 are repeated until a termination criterion is met. This criterion could be a maximum
number of iterations, reaching a specific function value threshold, or satisfying a desired precision.

Figure 4.8: Main steps of the downhill simplex algorithm in a 2-dimensional space. The
coefficients used in this example are not fixed and can vary.

The downhill simplex algorithm continues to iterate through these steps, gradually converging to-
wards the minimum (or maximum) of the function.

Limitations: It is a simple and robust optimization method, but it may converge slowly in some
cases or get trapped in local minima. This algorithm is advised up to 5 variables.

References
Nelder, John A, and Roger Mead. 1965. "A Simplex Method for Function Minimization." The Computer
Journal 7 (4): 308–13.

[Link]. Hooke-Jeeves Pattern Search

Available in: ModelCenter

Keywords— Mathematical programming

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 93
Optimization methods

Hooke-Jeeves Pattern Search is a non-gradient based algorithm for unconstrained [Link]

uses a simple search pattern where the current design is perturbed (plus and minus) in each dimen-
sion one at a time. If a better design is found, the current design point is moved to the new design.
Otherwise, the perturbation step size is reduced and the pattern search is repeated until the step
size becomes smaller than the convergence criteria.

[Link]. DAKOTA OPT++ Parallel Direct Search (PDS)

Available in: ModelCenter

Keywords— Mathematical programming

OPT++ PDS is an unconstrained optimization algorithm based on the Nelder-Mead simplex algorithm.

[Link]. DAKOTA Asynchronous Parallel Pattern Search (APPS)

Available in: ModelCenter

Keywords— Mathematical programming

The Asynchronous Parallel Pattern Search (APPS) (Gray and Kolda 2006) is a non-gradient based
optimization algorithm, which is a variation of the Hooke-Jeeves Pattern Search. In its non-blocking
(asynchronous) mode, pattern search moves to the next iteration as soon as a new trial point im-
proves the current design. In its blocking mode, all designs in the search pattern is evaluated before
moving to the next iteration. APPS uses a penalty function to handle nonlinear constraints. This
algorithm supports an asynchronous pattern search technique where the search along each offset
direction continues without waiting for searches along other directions to finish. The algorithm is
a derivative-free optimization method used to solve single-objective optimization problems. The
variable bounds are directly considered in search direction, while general nonlinear constraints are
handled by a penalty function. As the algorithm moves along search directions by a step length,
the step length is contracted in case of unsuccessful iterate. Constraint penalty and constraint tol-
erance factors can be specified to adjust constraint handling. The algorithm terminates when con-
vergence tolerance is reached or maximum function evaluations or maximum iterations is reached.

References
Gray, Genetha A, and Tamara G Kolda. 2006. "Algorithm 856: APPSPACK 4.0: Asynchronous Parallel
Pattern Search for Derivative-Free Optimization." ACM Transactions on Mathematical Software (TOMS)
32 (3): 485–507.

[Link]. DAKOTA Coliny DIRECT

Available in: ModelCenter

Keywords— Mathematical programming

The DIviding RECTangles (DIRECT) is a derivative free global optimization method for solving single-
objective optimization problems. This algorithm seeks to balance local search and global search. It
performs local search in promising regions and perform global search in under-explored regions
at the same time. DIRECT uses a fixed penalty factor for constraint violations. Example: As shown
in Figure 4.9 (p. 95), DIRECT adaptively subdivides the space of feasible design points so as to

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
94 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

guarantee that iterates are generated in the neighborhood of a global minimum in finitely many
iterations.

Figure 4.9: Subdivision steps of DIRECT algorithm (DAKOTA 5.0 Reference Manual, p. 70)

Algorithm selects points for analysis around already evaluated point by dividing design space into
smaller rectangles. The selected points once analyzed are compared and points around the best
are selected for evaluation. Depending on the results of the evaluation the algorithm selects points
in neighborhood of the new best point. This ensures the algorithm is pursuing the local minima.
Within the same iteration the algorithm selects points for iteration around the next best point from
previous evaluation. This ensures that the algorithm gets the global minimum. The algorithm con-
tinues the process and thus continuously reduces the design space to get in the neighborhood of
potential global minima.

DIRECT algorithm terminates with the standard MaxFunctionEvaluations and SolutionTarget specific-
ations. Additionally, the MaxBoxsizeLimit specification terminates the algorithm if the size of the
largest sub-region falls below this threshold. The MinBoxsizeLimit specification terminates DIRECT
algorithm if the size of the smallest sub-region falls below this threshold. In practice, this latter
specification is likely to be more effective at limiting DIRECT.s search.

[Link]. DAKOTA NCSU DIRECT

Available in: ModelCenter

Keywords— Mathematical programming

This is an unconstrained optimization algorithm similar to DAKOTA Coliny DIRECT and implemented
by the North Carolina State University (NCSU) (Gablonsky 2001). This is a derivative free optimization
method which works by dividing the design space into small rectangular sub-regions. This algorithm

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 95
Optimization methods

tries to balance between local search and global search to identify the best design for a multi-
modal objective function.

References
Gablonsky, J. 2001. "DIRECT Version 2.0 Userguide Technical Report CRSC-TR01-08." Center for Research
in Scientific Computation, North Carolina State University, Raleigh, NC.

[Link]. DAKOTA Coliny Pattern Search

Available in: ModelCenter

Keywords— Mathematical programming

Coliny Pattern Search is an extension of the traditional Hooke-Jeeves Pattern Search. It uses a set
of offsets as a search pattern and has several options to control the search patterns. This is a non-
gradient based optimization that may be suited for noisy functions. A simple penalty function is
used to handle constraint violations. Pattern search technique is a nongradient-based optimization
methods which use a set of offsets from the current iterate to locate improved points in the design
space. The Coliny pattern search technique includes a variety of specification components. Tradi-
tional pattern search methods search with a fixed pattern of search directions to try to find improve-
ments to the current iterate. The Coliny Pattern Search methods generalize this simple algorithmic
strategy to enable control of how the search pattern is adapted, as well as how each search pattern
is evaluated. Additional search pattern includes simplex with N+1 points.

4.2.4. Response surface based methods

Response surface methods replace the model responses by mathematical surrogate functions. On
this base, the optimization problem is solved using the approximation model instead of time consuming
solver calls. The classical procedure uses a Design of Experiments scheme which is evaluated by the
solver. Based on these designs often polynomial functions are used for the approximation (Myers and
Montgomery 2002).

Limitations: The polynomial degree does often not represent the nonlinearity of the solver model
and a low quality of the approximation leads often to unusable design suggestions. For this reason
it is not recommended to use global polynomial response surface methods only to search for an op-
timal design.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
96 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

Figure 4.10: Global polynomial approximation of the maximum amplitude (left, ) and
the damped eigen-frequency (right, ) of the optimized oscillator using a quadratic
basis and a full factorial design scheme

In order to demonstrate the weakness of global polynomial response surface approximation in com-
bination with classical DoE schemes, the maximum amplitude and the damped eigen-frequency of
the damped oscillator are approximated by a quadratic polynomial. Support points are generated
using a three-level full factorial design. Although the approximation quality is indicated to be very
well by means of the adjusted Coefficient of Determination, the optimum found on the approximation
( , ) is far away from the true optimum.

References
Myers, R., and D. C. Montgomery. 2002. Response Surface Methodology. 2nd ed. John Wiley & Sons,
Inc.

[Link]. Adaptive Response Surface Method (ARSM)

Available in: optiSLang

Keywords— Metamodel-based, Optimal Design

In order to improve the approximation quality around the optimum, adaptive methods are very
efficient. optiSLang provides a polynomial based local Adaptive Response Surface Method (ARSM).
The ARSM procedure starts at a single start design, and an initial Design of Experiments (DoE)
scheme is built having the start design as center point. For linear or quadratic polynomial models
the corresponding D-optimal design schemes are preferred as DoE schemes. Based on the approx-
imation of the model responses the optimal design is searched within the parameter bounds of
the DoE scheme. In the next iteration step a new DoE scheme is built around this optimal design.
Depending on the distance between the optimal designs of the current and previous iteration
steps, the DoE scheme is moved, shrunken or expanded. This procedure is shown in principle in
Figure 4.11 (p. 98).

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 97
Optimization methods

Further details about the adaptation procedure can be found in (Etman et al. 1996) and (Stander
and Graig 2002).

Figure 4.11: Adaptation of the polynomial approximation scheme inside the Adaptive Response
Surface Method

The ARSM algorithm converges if the DoE is shrunken to a minimum size or if the change of the
optimal design position and its objective value between two iteration steps is below a specified
tolerance. Since the DoE scheme uses 50% more designs than needed for the polynomial approx-
imation, the solver noise is smoothed and few failed designs are not problematic for the optimizer.
Due to its efficiency for up to 20 variables and its robustness against solver noise, the ARSM is the
method of choice for low dimensional single-objective optimization problems. Limitations: Applicable
up to 20 variables. It is recommended to use the optimal design of a preceding sensitivity analysis
as start design for the ARSM. For a strongly localized search the start range of the initial DoE scheme
should be reduced.

Example: In Figure 4.12 (p. 99) the convergence of the ARSM optimizer is shown for the oscillator
problem by analyzing the noisy objective function.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
98 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

Figure 4.12: Convergence of the Adaptive Response Surface Method for the damped oscillator
by using the noisy objective function: designs used for the adaptation (left) and modification
of the local DoE bounds during the iteration (right)

By using a start range of of the design space and a linear polynomial basis, the optimizer runs
within a few iteration steps in the region of the true optimum. After 20 iteration steps the algorithm
converges at , , where the constraint condition is fulfilled. This examples
clarifies that the ARSM optimizer shows a stable convergence behavior for noisy model responses
in contrast to the gradient based methods.

References
Etman, L., J. Adriaens, M. van Slagmaat, and A. Schoofs. 1996. "Crashworthiness Design Optimization
Using Multipoint Sequential Linear Programming." Structural Optimization 12: 222–28.

Stander, N., and K. Graig. 2002. "On the Robustness of a Simple Domain Reduction Scheme for
Simulation-Based Optimization." Engineering Computations 19: 431–50.

[Link]. Optimization using the Metamodel of Optimal Prognosis

Available in: optiSLang

Due to the power of the Metamodel of Optimal Prognosis (Section 3.4 (p. 68)) procedure in finding
an optimal variable subspace and an optimal approximation model for each investigated model
response, it is strongly recommended to use the MOP approximation for a first optimization step
instead of global polynomial models.

The Coefficient of Prognosis (Section 3.2.5 (p. 61)) gives a much more reliable estimate of the ap-
proximation quality than the Coefficient of Determination (Section 3.2.4 (p. 59)) as explained in
section Section 3.2.5 (p. 61). If the quality of a certain model response, which is used in the objective
or constraint functions, is low, the optimization procedure may not find a useful optimized design.
In order to check the quality of the optimal design found on the MOP, optiSLang provides a verific-
ation of the best design where the solver outputs and the true objective and constraint values are
calculated by a single solver call. Often the best design found by the MOP based optimization is

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 99
Optimization methods

better than the best design of the preceding sensitivity analysis. If this is the case, this design could
be used as a start design for a further local search.

Remark 1. Since the solver noise is smoothed by the MOP approximation, this procedure is more stable
than the gradient-based methods. Furthermore, non-convex optimization problems can be solved by
using global optimizers such as evolutionary algorithms (EA) on the approximation function.

Limitations: In the presence of many constraint conditions often the MOP based best design violates
one or more of these conditions if the approximation quality is not perfect. This problem can be
solved by adjusting the limits of the corresponding constraint conditions in order to push the op-
timizer running on the approximation back to the feasible region. After each adjustment of the
constraint conditions the best design should be verified to check for a possible fulfillment of the
critical constraints. Generally the MOP uses uniformly distributed designs for the approximation
model, which may lead to a insufficient local approximation quality around the optimum. Therefore,
optimization using global approximation models should be understood as a low-cost pre-optimiz-
ation step. Example: In order to demonstrate the benefit of an optimization on the MOP approx-
imation, the damped oscillator is investigated by building the MOP with 100 Latin Hypercube
samples. The approximation functions of the MOP are shown in Figure 4.13 (p. 100).

Figure 4.13: Approximation of the maximum amplitude (left, ) and of the damped
eigen-frequency (right, ) by using the Metamodel of Optimal Prognosis with 100
Latin Hypercube samples

The figure indicates an excellent approximation quality in terms of the Coefficient of Prognosis. By
using the MOP approximation functions the optimum is determined very close to the true optimum
( , ). The approximated damped eigen-frequency at the obtained optimum
is , but the real eigen-frequency verified by the solver is , which means that
the constraint condition is slightly violated. In this case a reduction of the maximum allowed eigen-
frequency would force the optimizer to stay in the feasible region.

[Link]. Adaptive Metamodel of Optimal Prognosis (AMOP)

Available in: optiSLang

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
100 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

Keywords— Metamodel-based, adaptive, MOP

By using the criteria-refinement of the AMOP algorithm a single and multi-objective optimization
task can by evaluated. Based on the estimated local approximation errors introduced in section
Section 3.4.2 (p. 71), the root mean squared prediction errors of the optimization criteria are estim-
ated for a large set of candidate points in the design space. Using the criterion of expected improve-
ment as in the EGO approach in section Section [Link] (p. 104) each candidate point can be evaluated
for a single-objective function as follows (Jones, Schonlau, and Welch 1998 (p. 101))

(4.10)

where and are the density and the distribution function of a unit Gaussian random variable.
The difference between the best solution found in the search so far, , and the approximated
objective function is standardized by the estimated approximation error using the
local estimate of the RMSE in Equation 3.93 (p. 72).

In case of a constraint optimization problem, the expected improvement of the objective function
in Equation 4.10 (p. 101) is multiplied by the probability of each inequality constraint condition
, that the candidate point is valid. By assuming again a normally distributed approximation
of the constraint functions, the constraint expected improvement reads

(4.11)

Limitations: Due to the internal use of Kriging and MLS approximation, the AMOP algorithm for
single-objective optimization is not advised for more than 10 input parameters.

References
Jones, Donald R, Matthias Schonlau, and William J Welch. 1998. "Efficient Global Optimization of
Expensive Black-Box Functions." Journal of Global Optimization 13 (4): 455.

[Link]. Adaptive Single-Objective

Available in: DesignXplorer, optiSLang

Keywords— LHS, Metamodel-based, Kriging, SQP

Adaptive Single-Objective (ASO) combines a latin hypecube sampling optimized (LHS optimized,
see Section [Link] (p. 14), also named Optimal Space Filling: OSF), a Kriging response surface
(Section 3.1.4 (p. 26)), and local optimizer (MISQP, see Section [Link] (p. 89)). It is a gradient-based
algorithm based on a response surface, which provides a refined, global, optimized result. ASO
supports a single objective and multiple constraints. It is available for continuous parameters, in-
cluding those with manufacturable values (list of continuous values). It does not support the use
of parameter relationships in the optimization domain and is available only for a direct approach.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 101
Optimization methods

Figure 4.14: ASO workflow

ASO starts with a first population (OSF), and refine and reduce the domain intelligently and auto-
matically in the next iterations to provide the global extrema.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
102 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

The Optimal Space-Filling Design is used for build a Kriging response surface. In the original OSF,
the number of samples equals the number of divisions per axis and there is one sample in each
division.

When a new OSF is generated after a domain reduction, the reduced OSF has the same number of
divisions as the original and keeps the existing design points within the new bounds. New design
points are added until there is a point in each division of the reduced domain.

In the following example (Figure 4.15 (p. 103)), the original domain has eight divisions per axis and
contains eight design points. The reduced domain also has eight divisions per axis and includes
two of the original design points. To have a design point in each division, six new design points
need to be added.

Figure 4.15: From the current OSF to the new OSF

Remark 1 (OSF of the reduced domain). The total number of design points in the reduced domain
can exceed the number in the original domain if multiple existing points wind up in the same division.
In the previous example above, if two existing points wound up in the same division of the new domain,
seven new design points (rather than six) would have been added to have a point in each of the remaining
divisions.

A response surface is created for each output, based on the current OSF and consequently on the
current domain bounds.

MISQP is run on the current Kriging response surface to find potential candidates. A few MISQP
processes are run at the same time, beginning with different starting points, and consequently,
giving different candidates.

All the obtained candidates are either validated or not, based on the Kriging error predictor. The
candidate point is checked to see if further refinement of the Kriging surface changes the selection
of this point. A candidate is considered as acceptable if there aren"t any points, according to this
error prediction, that call it into question. If the quality of the candidate is not called into question,
the domain bounds are reduced. Otherwise, the candidate is calculated as a verification point.

When a new verification point is calculated, it is inserted in the current Kriging response surface as
a refinement point and the MISQP process is restarted.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 103
Optimization methods

When candidates are validated, new domain bounds must be calculated. If all of the candidates
are in the same zone, the bounds are reduced, centered on the candidates. Otherwise, the bounds
are reduced as an inclusive box of all candidates. At each domain reduction, a new OSF is generated
(conserving design points between the new bounds) and a new Kriging response surface is generated
based on this new OSF. Limitations: Due to the internal use of Kriging, ASO is not advised for more
than 10 input parameters.

[Link]. Efficient Global Optimization

Available in: LS-OPT

Keywords— Metamodel-based, Bayesian Optimization, Expected Improvement

Bayesian global optimization (BO) techniques have been successfully used in various problems
(Mo kus 1975), (Mo kus 2005), (Jones, Schonlau, and Welch 1998). One of the most popular algorithms
is the so-called Efficient Global Optimization (EGO) algorithm of (Jones, Schonlau, and Welch 1998).
It consists in sampling the point that maximizes the so-called expected improvement (EI).

EGO is a technique used in global optimization problems to efficiently find the global optimum of
an objective function within a given search space. It is particularly useful when the objective function
is expensive to evaluate and may have multiple local optima.

The basic idea behind EGO is to iteratively build a surrogate model of the objective function using
a small number of function evaluations, and then use this surrogate model to guide the search for
the global optimum. The surrogate model is typically a statistical model, such as a Gaussian process,
that captures the behavior of the objective function based on the available data points. The EGO
algorithm proceeds in the following steps:

1. Initialization: a initial design of experiments (DoE), is sampled from the search space. The objective
function is evaluated at these points to obtain the corresponding function values.

2. Surrogate Model Construction: The DoE points and their corresponding function values are used
to train a surrogate model, which estimates the behavior of the objective function across the
search space. A Gaussian process is commonly used as the surrogate model due to its flexibility
and ability to model complex functions.

3. Acquisition Function: An acquisition function is defined to determine the next point to evaluate.
The acquisition function balances the exploration-exploitation trade-off by considering both
the uncertainty of the surrogate model predictions and the desirability of the function values.
Common acquisition functions include Expected Improvement (EI) and Probability of Improvement
(PI).

4. Selection of the Next Point: The acquisition function is maximized to select the next point to
evaluate. This point is chosen as the one that is expected to yield the highest improvement
over the current best solution.

5. Objective Function Evaluation: The objective function is evaluated at the selected point, and its
value is added to the available data.

6. Surrogate Model Update: The surrogate model is updated with the new data point to improve
its accuracy.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
104 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

7. Convergence Check: The algorithm checks if a convergence criterion is met. If not, it returns to
step 3 and continues the iteration. The convergence criterion can be based on the number of
iterations, the improvement in the objective function value, or other stopping conditions.

By iteratively updating the surrogate model and selecting points that are likely to improve the ob-
jective function, EGO effectively explores the search space and focuses the search towards promising
regions. This helps in efficiently finding the global optimum while minimizing the number of ex-
pensive evaluations of the objective function.

Let be a Gaussian process. Let further and denote respectively the mean and
the variance of the conditional process . Last, let be the minimum value of the response
on the sample where , that is . The EGO algorithm (Jones,
Schonlau, and Welch 1998) uses the expected improvement (Equation 4.12 (p. 105)) as sampling
criterion:

(4.12)

Using some Gaussian computations, can be explicitly computed:

(4.13)

Where and are respectively the cumulative and probability density functions of .
The EGO algorithm adds to the sample the point that maximizes .An illustration of 5 iterations
of EGO on a toy example is displayed in Figure 4.16 (p. 106).

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 105
Optimization methods

Figure 4.16: Example of 1-dimensional EGO algorithm. Dashed red line: real function, dashed
black line: value of the current minimum, solid line: Kriging prediction ( ), light blue area:
Kriging variance ( ), black squares: design points.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
106 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

References
Jones, Donald R, Matthias Schonlau, and William J Welch. 1998. "Efficient Global Optimization of
Expensive Black-Box Functions." Journal of Global Optimization 13 (4): 455.

Mo kus, Jonas. 1975. "On Bayesian Methods for Seeking the Extremum." In Optimization Techniques
IFIP Technical Conference: Novosibirsk, July 1–7, 1974, 400–404. Springer.

———. 2005. "The Bayesian Approach to Global Optimization." In System Modeling and Optimization:
Proceedings of the 10th IFIP Conference New York City, USA, August 31–September 4, 1981, 473–81.
Springer.

[Link]. Probabilistic Inference for Bayesian Optimization

Available in: optiSLang

Keywords— Metamodel-based, Bayesian Optimization, Expected Improvement

Probabilistic Inference for Bayesian Optimization (PI-BO) is an algorithm based on the so-called
Bayesian Optimization (Jones, Schonlau, and Welch 1998). This type of optimization exploits the
computable model uncertainty of probabilistic machine learning models such as DIM-GP (Sec-
tion 3.1.10 (p. 53)) to propose new designs that are most likely to improve the objective functions
and satisfy any constraints that may exist. This usually involves starting with a smaller number of
initial designs and building an initial model.

Example: See Figure 4.17 (p. 108), where two starting designs were chosen and an initial DIM-GP
model was trained (green line). The model uncertainty is marked in gray, here the confidence
interval of a Gaussian distribution. The black dashed line represents the function to be maximized.
Now, instead of optimizing the objective function directly, a so-called acquisition function is optim-
ized (right image). There are several acquisition functions that can be used in Bayesian Optimization.
Most of them always represent a combination of probably the greatest improvement of the objective
function and reduction of the model uncertainty. This leads to the fact that there is always an al-
ternation between improvement of the objective function and exploration of the design space. In
Figure 4.17 (p. 108). the course of the acquisition function is shown on the right. This is maximized

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 107
Optimization methods

at the right edge of the function where, considering the mean of the function (green line) and the
model uncertainty (gray area), there is the greatest potential to improve the objective function.

Figure 4.17: Beginning of Bayesian Optimization.

i.e., that would be the next design that is evaluated in the simulation or experiment and then made
available to the model again. The model performs a new training and adapts to the new data. This
changes the mean of the model and also the model uncertainty. This process is repeated until
either the optimum is found (convergence) or the design budget is used up. Figure 4.18 (p. 108)
shows how the model looks after 6 or 10 iterations of this optimization and where the individual
designs were proposed. After 10 iterations the global maximum of the objective function was found
with multiple designs around it.

Figure 4.18: Illustration of the Bayesian optimization iterations and the associated acquisition
function.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
108 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

This type of optimization is a very efficient way of adaptive design of experiments where a probab-
ilistic machine learning model like DIM-GP helps the user to choose the next experiments and learns
continuously based on the intermediate results. For failed designs or designs that violate constraints,
the model will learn to avoid these areas in the future. This saves unnecessary design evaluations
in non-permissible areas. Of course, it is also possible to propose more than one design per adapt-
ation, depending on the possibility of evaluating the designs in parallel, this can be useful.

References
Jones, Donald R, Matthias Schonlau, and William J Welch. 1998. "Efficient Global Optimization of
Expensive Black-Box Functions." Journal of Global Optimization 13 (4): 455.

4.2.5. Nature-inspired methods

Nature-inspired optimization methods imitate natural processes like biological evolution or swarm
intelligence. Based on the principle "survival of the fittest" a population of artificial individuals searches
the design space of possible solutions in order to find a better approximation for the solution of the
optimization problem. All kinds of nature-inspired algorithms are following a particular paradigm. It
starts with a random initialization (see Figure 4.19 (p. 110) step 1) of a population containing individuals
representing possible solutions to the given optimization problem. After initialization the actual iter-
ation loop starts, where every loop represents a generation. To create a new generation the fittest
designs will be selected (step 2). The adaption referred to step 3 can be done in different ways, e.g.
crossover or swarm movement. Mutation (step 4) can be applied afterwards. All individuals are eval-
uated (step 5) by assigning them a fitness value based on the objective value and possible constraint
violations. Thus a high fitness score means a good accommodation to the problem and corresponds
with a small objective value. After evaluating each design the archive, which keeps good solutions,
will be updated for the next selection step (step 6). The generation counter is increased and the iter-
ation continues until a stopping criterion is fulfilled (step 7).

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 109
Optimization methods

Figure 4.19: General flowchart of nature-inspired optimization algorithms

The efficiency of nature-inspired methods can be significantly improved by choosing a suitable start
population. If the optimization is performed after an initial sensitivity analysis, it is strongly recom-
mended to use the best designs of the sensitivity analysis as start population. In order to keep the
global character of the optimization it is useful to select designs covering a certain range of the design
space. .

The usage of nature-inspired algorithms is recommended in a wide range of applications. They are
often the last resort for solving mathematically ill conditioned problems since their performance may
be still good in case of a high number of design variables or a high amount of failed designs. Never-
theless, the convergence behavior may be slower compared to other optimization methods. The use
of nature-inspired algorithms is recommended in all cases where gradient based optimization or re-
sponse surface approximation fails, in case of a high number of variables or constraints, in case of
discrete or binary design variables, in case of discrete responses or if the user is unaware about the
optimization problem.

[Link]. Evolutionary algorithms

Available in: ModelCenter, optiSLang

Keywords— Natures-inspired, Evolutionary Algorithm

Evolutionary algorithms (EA) are stochastic search methods that mimic processes of natural biolo-
gical evolution. Many EA variants have been implemented over the past decades, based on the
evolution strategies (ES) (Rechenberg 1964) and evolutionary programming (EP) (Fogel, Owens, and
Walsh 1966). These algorithms have been originally developed to solve optimization problems
where no gradient information is available, like binary or discrete search spaces, although they can
also be applied to problems with continuous variables. Within optiSLang there is a flexible imple-
mentation of genetic algorithms and evolutionary strategies available which allows to scale between
pure (strong) GA and pure (strong) ES.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
110 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

After the random or manual initialization of a start population the actual iteration loop starts, where
every loop represents a generation . First the best individuals based on their fitness values are
determined for reproduction. This parent selection sets the direction of the stochastic search process.
Out of this pool for reproduction individuals are selected by stochastic selection operators. These
operators are provided in order to adjust the selection pressure. Variation is introduced by applying
crossover and mutation operators to the selected individuals. These new individuals are called
newborns or offspring. Finally the resulting offspring individuals are evaluated and a new loop
starts.

The main difference between the two main variants, genetic algorithms (GA) and evolution strategies
(ES), is the way variation is introduced to the population. Recombination of genes using crossover
operators represents the main variation within genetic algorithms, while for evolution strategies
(adaptive) mutation introduces variation to the population (see Figure 4.20 (p. 111)).

Figure 4.20: Evolution strategies versus genetic algorithms

The selection of individuals for reproduction is based on random selection methods and requires
an assignment of fitness values. A rank-based fitness assignment is used to overcome the scaling
problem of the proportional fitness assignment. Due to the random selection and the fitness ranking,
the probability of selecting a design with low fitness is low. However, this is not impossible, so also
a weaker design could be used for the generation of offspring.

The crossover operator is a method of recombination where two parent individuals produce two
offspring by sharing information between chromosomes. The intention is to obtain individuals with
better characteristics (exploitation) and to maintain the diversity of the population (exploration).
Crossover is regarded as the main search operator in genetic algorithms.

Mutation introduces random variation to the genes of the offspring chromosome. Each gene is
selected for mutation with a specified probability or mutation rate respectively. The real-valued

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 111
Optimization methods

mutation is based on a normal distribution function for each gene with the value of the gene as
its mean value. Mutation is the main search operator for evolutionary strategies (ES) but can also
be applied after recombination (see Figure 4.19 (p. 110)). The mutation rate and the standard deviation
defined in relation to the variable range are the control parameters of the mutation which can stay
constant or get modified during the run of the algorithm. optiSLang provides constant and adaptive
mutation procedures. Further details can be found in (Bäck 1996) and (Riedel et al. 2005).

The final operator, the archive update scheme, specifies how the population of the next generation
is formed out of the best individuals of the old population and the generated offspring individuals.
The best individuals for the next generation are determined from a pool of individuals. Different
strategies of constituting this pool are available by adjusting the size of the archive ( ).

• ( )-strategy: archive size = 0

the best out of offspring form the parents for the next generation

• ( )-strategy:

the best individuals out of the best from the old population plus the offspring form the
parents for the next generation

• ( )-strategy: archive size =

the best individuals are determined out of the best from the old population plus the off-
spring

The ( )-strategy with empty archive leads to a complete replacement of the population at every
generation step with a maximum life span of an individual of only one generation. This method is
used for classic genetic algorithms. The ( )-strategy prevents parent individuals from being re-
placed by offspring with worse fitness. This strategy is recommended for use with evolution
strategies. The ( )-strategy allows the adjustment of the archive update scheme between
the two extreme cases.

optiSLang provides two constraint handling methods taking into account all individuals of the
population, which are compared regarding their objective values and constraint violations. The in-
feasible solutions are either penalized according to the scaled and summed values of violation or
to the number of constraint violations. Their fitness is determined by adding the fitness of the worst
feasible solution to the penalty value. These methods can handle multiple constraints and even if
all individuals are violating constraints, they can guide the search towards a feasible region of the
constraint space.

In optiSLang a predefined global and local search is available for the evolutionary algorithms. The
global search starts with a random or manually selected start population and is suitable to detect
new possible solutions in the whole design space. The local search tries to improve a single design
without discovering new regions.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
112 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

Figure 4.21: Convergence of the evolutionary algorithm with global search (left) and local
search (right) for the damped oscillator with noisy objective function

In Figure 4.21 (p. 113) both search strategies are compared. The figure indicates, that the global
search, which starts with a random start population, generates designs in the whole design space
while the designs generated with the local search method stay close to the start design.

References
Bäck, T. 1996. "Evolution Strategies: An Alternative Evolutionary Algorithm." Lecture Notes in Computer
Science 1063/1996: 1–20.

Fogel, L. J., A. J. Owens, and M. J. Walsh. 1966. Artificial Intelligence Through Simulated Evolution.
John Wiley & Sons.

Rechenberg, I. 1964. "Kybernetische Lösungsansteuerung Einer Experimentellen Forschungsaufgabe."

Riedel, J., S. Blum, R. Puisa, and M. Wintermantel. 2005. "Adaptive Mutation Strategies for Evolutionary
Algorithms: A Comparative Benchmark Study." In Proc. Weimarer Optimierungs- Und Stochastiktage
2.0, Weimar, Germany.

[Link]. Darwin algorithm

Available in: ModelCenter

Keywords— Natures-inspired, Evolutionary Algorithm

Darwin is a genetic search algorithm developed specifically for solving engineering optimization
problems. Darwin is capable of handling discrete variables, continuous variables, and any number
of constraints. Because Darwin does not require gradient information, it is able to effectively search
non-linear and noisy design spaces.

Darwin is a genetic optimizer that can solve constrained design problems. It supports continuous
and discrete variables. A penalty function is used to handle violated constraints. It uses elitist
method to keep best design(s) from the previous generation.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 113
Optimization methods

[Link]. EVOLVE
Available in: ModelCenter

Keywords— Natures-inspired, Evolutionary Algorithm

EVOLVE is a genetic search algorithm that can be used to optimize problems with a mix of continu-
ous, integer and discrete design variables. Constraints are considered by using an exterior penalty
function. It has a capability to improve global search (called sharing) by penalizing designs that are
near to other good designs.

[Link]. DAKOTA Coliny Evolutionary Algorithm (EA)

Available in: ModelCenter

Keywords— Natures-inspired, Evolutionary Algorithm

Coliny Evolutionary Algorithm is a very flexible implementation of the genetic algorithm. This is a
derivative free algorithm that handles both discrete and continuous design variables. It uses a
simple penalty scheme to handle constraint violations. Basic steps of evolutionary algorithm are as
follows:

1. Select an initial population randomly and perform function evaluations on these individuals

2. Perform selection for parents based on relative fitness

3. Apply crossover and mutation to generate new individuals from the selected parents

• Apply crossover with a fixed probability from two selected parents

• If crossover is applied, apply mutation to the newly generated individual with a fixed probab-
ility

• If crossover is not applied, apply mutation with a fixed probability to a single selected parent

4. Perform function evaluations on the new individuals

5. Perform replacement to determine the new population

6. Return to step 2 and continue the algorithm until convergence criteria are satisfied or iteration
limits are exceeded

The Evolutionary algorithm decides on the initial population based on the InitializationType. The
algorithm generates up to PopulationSize number of designs randomly at the start. Once population
is fixed, the FitnessType option guides the process of selection of parents for crossover. The Crossov-
erType controls what approach is employed for combining parent genetic information to create
offspring, and the CrossoverRate specifies the probability of a crossover operation being performed
to generate a new offspring. The MutationType option controls what approach is employed in ran-
domly modifying continuous design variables within the Evolutionary Algorithm population. Once
crossover is complete for the parent designs the ReplacementType option helps the algorithm decide
how current population and newly generated solutions would combine to form the new population.
Thus the algorithm completes one cycle. Algorithm again performs selection based on relative fitness
to loop through the steps until one of the termination criterions is encountered.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
114 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

[Link]. Covariance Matrix Adaptation

Available in: ModelCenter, optiSLang

Keywords— Natures-inspired, Evolution Strategy Covariance Matrix Adaptation (CMA) is a

stochastic, or randomized, evolution strategy for real-parameter (continuous domain) optimization
of non-linear, non-convex functions. It is designed to tackle ill-conditioned and non-separable
problems. The start population is either given as start designs or randomly generated by a uniformly
distributed Latin Hypercube Sampling (Section 2.2.2 (p. 12)) within the defined input parameter
ranges and may be smaller or larger than the actual population size. After evaluating the start
population, the best designs are determined. They are used to find the mean value for the
first sampling of the offsprings. From that point on in each iteration the best designs are de-
termined and used to update the mean value and the covariance matrix of the distribution new
candidate solutions are sampled according to. In contrast to most classical methods, fewer assump-
tions on the nature of the underlying objective function are made. Only the ranking between can-
didate solutions is exploited for learning the sample distribution and neither derivatives nor even
the function values themselves are required by the method. Limitations: This method is advised
up to 20 variables. Example: In Figure 4.22 (p. 115) the convergence of the CMA approach is shown
for the optimization of the Rosenbrock function in 2D. Every group of colored points represents a
population. The true optimum (1,1) is achieved for this example.

Figure 4.22: Convergence of the covariance matrix adaptation for the rosenbrock function

[Link]. Particle Swarm Optimization

Available in: ModelCenter, optiSLang

Keywords— Natures-inspired, Particle Swarm

Particle Swarm Optimization (PSO) is a nature-inspired optimization method developed by (Kennedy

and Eberhart 1995) that imitates the social behavior of a swarm. Information about preferred positions
will be transmitted to other individuals of the swarm such that they move into directions of previous
optimal positions. In nature an individual can be represented by a bird, bee or fish. An overview
of modifications and improvements of this class of algorithms can be found e. g. in (Engelbrecht
2005) and (Poli, Kennedy, and Blackwell 2007).

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 115
Optimization methods

Figure 4.23: Update of a particle position in the Particle Swarm Optimization using a
combination of the old velocity , the direction to the local best position and the direction
to the global best position

PSO starts with a random or manual initialization of a start population which may be larger than
the actual population size . Out of this start population the best individuals are determined
where each one represents a possible solution of the optimization problem. Swarm intelligence is
influenced by two main components representing the personal and global behavior. This means
that each individual remembers its personal best solution and will also be influenced by the best
solution of the swarm. Each individual will change its position into the direction of its personal
best found position and the global best found position :

(4.14)

There are three important parameters that influence the speed and spread of the swarm. The inertia
weight is a scaling factor for the velocity of the previous iteration step . The personal acceleration
coefficient is a scaling factor for the second term, which is also called cognitive component. The
third term, called social component, is scaled by the global acceleration coefficient . and
are two vectors of random numbers uniformly chosen from .

The choice of suitable coefficients , and is very important for the convergence behavior of
the PSO. optiSLang provides two predefined search strategies with default values for each coefficient.
A local search strategy is recommended if the user has preliminary information about the optimiz-
ation space and has already found a pre-optimized design. In this case the swarm movement is
slower and less intensive throughout all generations. In a predefined global search strategy the
swarm movement is very intensive at the beginning and will be damped throughout the optimization
process by decreasing the weights of the previous velocity and the cognitive component and in-

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
116 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

creasing the weight of the social component. In that way a global exploration is reached in the first
iterations and exploitation is possible in last iterations.

The global best solution is recorded in the archive and is updated at the end of each iteration step.
In the presence of constraint conditions the global best solution is found by taking the feasible
solution with minimal objective value. If there are only infeasible solutions in the first population
the global best will taken as the individual with least constraint violations.

The PSO works similarly to the evolutionary algorithms for optimization problems with a large
amount of failed designs. For optimization tasks with continuous design variables the PSO converges
generally very close to the global optimum. Similar to evolutionary algorithms, the choice of a
qualified start population e. g. from sensitivity analysis may significantly improve the efficiency of
the algorithm. Limitations: PSO is less efficient compared to evolutionary algorithms if discrete
design variables and many constraint conditions are investigated.

References
Engelbrecht, A. P. 2005. Fundamentals of Computational Swarm Intelligence. John Wiley & Sons.

Kennedy, J., and R. Eberhart. 1995. "Particle swarm optimization."

Poli, R., J. Kennedy, and T. Blackwell. 2007. "Particle Swarm Optimization, An Overview." Swarm In-
telligence 1: 33–57.

[Link]. Stochastic Design Improvement

Available in: optiSLang

Keywords— Natures-inspired, Stochastic Approach

Stochastic Design Improvement (SDI) is a local, single-objective optimization procedure that improves
a proposed design by using a simple stochastic approach without having extensive knowledge
about potential relationships in design space. Based on a start population which is either given as
start designs or randomly generated by a uniformly distributed Latin Hypercube Sampling within
the defined input parameter ranges, the first best design is determined. In each iteration step, the
best design is taken as new center point for the sampling of the next iteration as shown in Figure
4.24 (p. 118). The ranges of the sampling scheme are adapted in each iteration step. Depending on
the optimization problem the whole population might move into a better region and achieve an
improvement in each step. All individuals of one iteration will be compared concerning objective
and constraints by using the parameter-free constraint handling method applied in the evolutionary
algorithms. The algorithm converges if either a maximum number of iterations was reached or if
for a specified number of iterations there was no improvement.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 117
Optimization methods

Figure 4.24: Update scheme of the Stochastic Design Improvement approach by generating
a random sampling around the best design of a previous iteration

Due to its pure stochastic approach, the SDI is very robust and works for a high amount of failed
designs, discrete and continuous parameters, and a high number of design variables and constraint
conditions.

Limitations: SDI was designed to improve an initial design but not to find global optimal solutions.
For this reason its efficiency may be significantly lower compared to evolutionary algorithms and
Particle Swarm Optimization.

[Link]. Adaptive Simulated Annealing

Available in: LS-OPT

Keywords— Natures-inspired, Simulated Annealing

The Simulated Annealing (SA) is a global stochastic optimization algorithm that mimics the metal-
lurgical annealing process. The Simulated Annealing process starts with an initial solution and then
iteratively improves the current solution by randomly perturbing it and accepting the perturbation
with a certain probability. The probability of accepting a worse solution is initially high and gradually
decreases as the number of iterations increases.

The original simulated annealing algorithm was developed as a generalization of the Metropolis
Monte Carlo integration algorithm (Metropolis et al. 1953) to solve various combinatorial problems
by Kirkpatrick et al. (Kirkpatrick, Gelatt Jr, and Vecchi 1983). The term "simulated annealing" derives
from the rough analogy of the way that the liquefied metals at a high temperature crystallize on
freezing. At high temperatures, the atoms in the liquid are at a high energy state and move freely.
When the liquid is cooled, the energy of the molecules gradually reduces as they go through many
lower energy states, and consequently their motion. If the liquid metal is cooled too quickly or
"quenched", the atoms do not get sufficient time to reach thermal equilibrium at a temperature
and might result in a polycrystalline structure with higher energy. This atomic structure of material
is not necessarily the most desired. However, if the rate of cooling is sufficiently slow, the atoms

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
118 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

are often able to achieve the state of minimum (most stable) energy at each temperature state,
resulting in a pure crystalline form. This process is termed as "annealing". Kirkpatrick et al.(Kirkpatrick,
Gelatt Jr, and Vecchi 1983) employed this annealing analogy to develop an efficient search algorithm.
Pincus (Pincus 1970), and Cerny ( ern 1985) also are also independently credited with the develop-
ment of modern simulated annealing algorithm. In simulated annealing parlance, the objective
function of the optimization algorithm is often called "energy" and is assumed to be related to
the state, popularly known as temperature , by a probability distribution. The Boltzmann distribution
is the most commonly used probability distribution:

Probability(E) ,

where is the Boltzmann"s constant.

[Link].1. Algorithm
The search initializes with the temperature being high and cooling slowly such that the system
goes through different energy states in search of the lowest energy state that is the global minima
of the optimization problem. A stepwise description of the simulated annealing algorithm is as
follows:

1. Initialization: The search starts by identifying the starting state and corresponding
energy . The temperature is initialized at a high value: . A cooling
schedule, acceptance function, and stopping criterion are defined. This is iteration
.

2. Sampling: A new point is sampled using the candidate distribution , and set
, and corresponding energy is calculated

3. Check acceptance: Sample a uniform random number in [0, 1] and set if

or otherwise. Where is the acceptance function that determ-
ines if the new state is accepted.

4. Temperature update: Apply the cooling schedule to the temperature: .

5. Convergence check: Stop the search if the stopping criterion is met, else set and go
to Step 2.

As is obvious, the efficiency of the simulated annealing algorithm depends on appropriate choices
of the mechanism to generate new candidate states D, cooling schedule C, acceptance criterion
A, and stopping criterion. While many options have been proposed in literature, the very fast
simulated reannealing methodology proposed by Ingber (1989) (Ingber 1989) has been the most
promising. This algorithm is also known as adaptive simulated annealing (ASA) (Ingber et al.
1993). The different selections along with a very brief historical perspective are outlined as follows.

[Link].2. Acceptance function

Two most prominent acceptance functions used to accept a candidate point are the Metropolis
criterion and the Barker criterion. Metropolis criterion : .

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 119
Optimization methods

Barker criterion : . The theoretical motivation for such a

restricted choice of acceptance functions can be found in (Schuur 1997). It is also shown that
under appropriate assumptions, many acceptance functions, which share some properties, are
equivalent to the above two criteria. The Metropolis criterion is the most commonly used selection
criterion and this is chosen as the acceptance function in LS-OPT.

[Link].3. Sampling algorithm

The choice of the next candidate distribution and the cooling schedule for the temperature are
typically the most important (and strongly interrelated) issues in the definition of a SA algorithm.
The next candidate state, , is usually selected randomly among all the neighbors of the current
solution, , with the same probability for all neighbors. The choice of the size of the neighborhood
typically follows the idea that when the current function value is far from the global minimum,
the algorithm should have more freedom, i.e., larger "step sizes" are allowed. However, Ingber
(Ingber 1989) suggested using a more complicated, non-uniformly selection procedure outlined
below to allow much faster cooling rates. Let design variable be bounded as, . Then
the new sample is given by:
(4.15)

Where is estimated as follows:

(4.16)

The most important distinction in ASA with standard SA is the use of an independent temperature
schedule ( ) for each parameter along with the temperature associated with the energy function.
The cooling schedule for the parameter temperature, used to generate N dimensional design
vector, is:
(4.17)

(4.18)

The control parameter depends on two free parameters and , defined as

(4.19)

The ratio is the parameter temperature ratio and the parameter is linked to
the time allowed (number of steps) at each parameter temperature state. Ingber (Ingber 2000)
found that the search procedure is sensitive to the choice of the two parameters and should be
selected carefully. Relatively, the parameter temperature ratio is the more important of the two
parameters.

[Link].4. Cooling schedule

The basic idea of the cooling schedule is to start the algorithm off at high temperature and then
gradually drop the temperature to zero. The primary goal is to quickly reach the so called effective
temperature, roughly defined as the temperature at which low function values are preferred but
it is still possible to explore different states of the optimized system (Basu and Frazer 1990). After
that the simulated annealing algorithm lowers the temperature by slow stages until the system
"freezes" and no further changes occur. Geman and Geman (Geman and Geman 1984) found the

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
120 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

lower bound on the cooling schedule to be where is an artificial time measure of the
annealing schedule. Hence,
(4.20)

This strategy is also known as Boltzmann annealing (Szu and Hartley 1987). Later this strategy
was modified (Laarhoven and Aarts 1987) to enable a much faster cooling schedule of:
(4.21)

A straightforward and most popular strategy is to decrement by a constant factor every it-
erations:
(4.22)

where is slightly greater than 1 (e.g. = 1.001). The value of should be large enough, so
that "thermal equilibrium" is achieved before reducing the temperature. A rule of thumb is to
take proportional to the size of neighborhood of the current solution. Nevertheless, the fastest
cooling rate was made possible by using Ingber"s algorithm that allowed an exponentially faster
cooling rate of
(4.23)
(4.24)

As was described in the previous section, the cooling rate is governed by the two free parameters
that are linked to the temperature ratio and annealing scale,
(4.25)

Typically the temperature ratio used to drive the energy (objective) function is linked to the
parameter temperature ratio called here as "cost-parameter annealing ratio".

[Link].5. Stopping criterion

Selection of an appropriate stopping criterion is one of the most difficult tasks in stochastic op-
timization algorithms because it is unknown a priori if the algorithm has reached the global optima
or is stuck in a hard local optimum. Thus the stopping rules proposed in the literature about SA,
all have a heuristic nature and are, in fact, more problem dependent than algorithm dependent.
Some common ideas in the heuristics are i) stop when it does not make a noticeable progress
over a number of iterations, ii) stop when the number of function evaluations is reached, and iii)
stop when the temperature has fallen substantially to a desired minimum level . The last two
criteria are used to terminate the adaptive simulated annealing search in LS-OPT.

[Link].6. Re-annealing
For multi-dimensional problems, most often the objective function has variable sensitivities for
different parameters and at different sampling states. Hence, it is worth while to adjust the
cooling rates for different parameters. Ingber (Ingber 1989) used a reannealing algorithm to
periodically update the annealing time associated with parameters and the energy function such
that the search is more focused in the regions with potential of improvements. For this, he sug-
gested computing the sensitivities of the energy function as,
(4.26)

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 121
Optimization methods

All the annealing time parameters are updated by the largest sensitivity as follows:
(4.27)

(4.28)

The new annealing time associated with the parameter is . Similarly the temperature
parameter associated with the energy function is scaled. One can easily deduce from the above
formulation that reannealing stretches the ranges of the insensitive parameters relative to the
sensitive parameters. More details of reannealing can be obtained elsewhere (Ingber 2000).

[Link].7. Some comments

1. It is difficult to find the initial temperature directly, because this value depends on the
neighborhood structure, the scale of the objective function, the initial solution, etc. In (Kirk-
patrick, Gelatt Jr, and Vecchi 1983) a suitable initial temperature is one that results in an average
uphill move acceptance probability of about 0.8. This can be estimated by conducting
an initial search, in which all uphill moves are accepted and calculating the average objective
increase observed. In some other papers, it is suggested that parameter is set to a value,
which is larger than the expected value of that is encountered from move to move.
In (Basu and Frazer 1990) it is suggested to spend most of the computational time in short
sample runs with different in order to detect the effective temperature. In practice, the
optimal control of may require physical insight and trial-and-error experiments. According
to (Bounds 1987), "choosing an annealing schedule for practical purposes is still something
of a black art".

2. Simulated annealing has proved surprisingly effective for a wide variety of hard optimization
problems in science and engineering. Many of the applications in our list of references attest
to the power of the method. This is not to imply that a serious implementation of simulated
annealing to a difficult real world problem will be easy. In the real-life conditions, the energy
trajectory, i.e. the sequence of energies following each move accepted, and the energy land-
scape itself can be highly complex. Note that state space, which consists of wide areas with
no energy change, and a few "deep, narrow valleys", or even worse, "golf-holes", is not suited
for simulated annealing, because in a "long, narrow valley" almost all random steps are uphill.
Choosing a proper stepping scheme is crucial for SA in these situations. However, experience
has shown that simulated annealing algorithms are more likely trapped in the largest basin,
which is also often the basin of attraction of the global minimum or of the deep local minimum.
Anyway, the possibility, which can always be employed with simulated annealing, is to adopt
a multi-start strategy, i.e. to perform many different runs of the SA algorithm with different
starting points.

3. Another potential drawback of using SA for hard optimization problems is that finding a good
solution can often take an unacceptably long time. While SA algorithms may quickly detect
the region of the global optimum, they often require a few iterations to improve its accuracy.
For small and moderate optimization problems, one may be able to construct effective pro-
cedures that provide similar results much more quickly, especially in cases when most of the
computing time is spent on calculations of values of the objective function. However, it should
be noted that for large-scale multidimensional problems an algorithm which always (or often)
obtains a solution near the global optimum is valuable, since various local deterministic op-
timization methods allow quick refinement of a nearly correct solution.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
122 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

In summary, simulated annealing is a powerful method for global optimization in challenging

real world problems. Certainly, some "trial and error" experimentation is required for an effective
implementation of the algorithm. The energy (cost) function should employ some heuristic related
to the problem at hand, clearly reflecting how "good" or "bad" a given solution is. Random per-
turbations of the system state and corresponding cost change calculations should be simple
enough, so that the SA algorithm can perform its iterations efficiently. The scalar parameters of
the simulated annealing algorithm have to be chosen carefully. If the parameters are chosen such
that the optimization evolves too fast, the solution converges directly to some, possibly good,
solution depending on the initial state of the problem.

References

Basu, Atanu, and L Neil Frazer. 1990. "Rapid Determination of the Critical Temperature in Simulated
Annealing Inversion." Science 249 (4975): 1409–12.

Bounds, David G. 1987. "New Optimization Methods from Physics and Biology." Nature 329 (6136):
215–19.

ern , Vladimı r. 1985. "Thermodynamical Approach to the Traveling Salesman Problem: An Efficient
Simulation Algorithm." Journal of Optimization Theory and Applications 45: 41–51.

Geman, Stuart, and Donald Geman. 1984. "Stochastic Relaxation, Gibbs Distributions, and the
Bayesian Restoration of Images." IEEE Transactions on Pattern Analysis and Machine Intelligence,
no. 6: 721–41.

Ingber, Lester. 1989. "Very Fast Simulated Re-Annealing." Mathematical and Computer Modelling
12 (8): 967–73.

Ingber, Lester et al. 1993. "Adaptive Simulated Annealing (ASA)." Global Optimization C-Code,
Caltech Alumni Association, Pasadena, CA.

Ingber, Lester. 2000. "Adaptive Simulated Annealing (ASA): Lessons Learned." arXiv Preprint
Cs/0001018.

Kirkpatrick, Scott, C Daniel Gelatt Jr, and Mario P Vecchi. 1983. "Optimization by Simulated An-
nealing." Science 220 (4598): 671–80.

Laarhoven, Peter J. M. van, and Emile H. L. Aarts. 1987. Simulated Annealing. Dordrecht: Springer
Netherlands.

Metropolis, Nicholas, Arianna W Rosenbluth, Marshall N Rosenbluth, Augusta H Teller, and Edward
Teller. 1953. "Equation of State Calculations by Fast Computing Machines." The Journal of Chemical
Physics 21 (6): 1087–92.

Pincus, Martin. 1970. "A Monte Carlo Method for the Approximate Solution of Certain Types of
Constrained Optimization Problems." Operations Research 18 (6): 1225–28.

Schuur, Peter C. 1997. "Classification of Acceptance Criteria for the Simulated Annealing Algorithm."
Mathematics of Operations Research 22 (2): 266–75.

Szu, Harold, and Ralph Hartley. 1987. "Fast Simulated Annealing." Physics Letters A 122 (3-4):
157–62.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 123
Optimization methods

[Link]. Differential Evolution

Available in: Model Center

Keywords— Natures-inspired, Direct Search

Brief description: Differential Evolution (DE) is a parallel direct search method developed by Storn
and Price in 1997 (Storn and Price 1997) and primarily targeted at unconstrained global optimization
problems. It is similar to the Genetic Algorithm and Particle Swarm Optimization in that it uses a
population of parameter vectors which evolves over time. DE follows the same basic steps also
employed generation-wise in the GA, namely Mutation, Crossover and Selection, although the detail
of the implementation is quite different.

The basic steps are described in Reference (Storn and Price 1997) as follows:

Mutation

For each target vector the population size, a mutant vector is generated ac-
cording to
(4.29)

With random indices , integer, mutually different and . The randomly

chosen integers are also chosen to be different from the running index , so that
must be to allow for this condition. is a real and constant factor [0,2] which controls the
amplification of the differential variation .

Crossover
To increase the diversity of the perturbed parameter vectors, crossover is introduced. To this end
the trial vector:
(4.30)

is formed where

(4.31)

In the last equation, is the th evaluation of a uniform random number generator with
outcome [0,1]. is the crossover constant [0,1] which has to be determined by the user.
is a randomly chosen index which ensures that gets at least one parameter
from .

Selection

To decide whether or not it should become a member of generation +1, the trial vector is
compared to the target vector using the greedy criterion. If vector yields a smaller cost
function value than the is set to ; otherwise the old value is retained.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
124 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

Other DE variants
In (Storn and Price 1997) Storn and Price also discuss other variants of DE in which they use the
notation: DE/x/y/z where

• x specifies the vector to be mutated which can be "rand" (a randomly chosen population vector)
or "best" (the vector of lowest cost from the current population).

• y is the number of difference vectors used.

• z denotes the crossover scheme. The variant described above is known as "bin" (crossover due
to independent binomial experiments).

Hence the strategy as described above is known as DE/rand/1/bin. The paper also mentions that
DE/best/2/bin could be beneficial. In this case
(4.32)

It is concluded that the usage of two difference vectors seems to improve the diversity of the
population if the number of population vectors is high enough. Because DE is such a simple
strategy, the astonishing result from (Storn and Price 1997) was that several examples showed
DE/rand/1/bin to be superior to several other algorithms tested, namely Adaptive Simulated Annealing
(Ingber et al. 1993), the Annealed Nelder and Mead approach (Press et al. 1992), the Breeder Genetic
Algorithm (Mühlenbein and Schlierkamp-Voosen 1993), the EASY Evolution Strategy (Voigt 2005),
and the method of Stochastic Differential Equations (Aluffi-Pentini, Parisi, and Zirilli 1985).

Constrained Optimization
The DE algorithm implemented in LS-OPT was developed by Kitayama et al. (Kitayama, Arakawa,
and Yamazaki 2011). The differences between the Kitayama and Storn & Price algorithms are mostly
of a minor nature. Whereas Storn and Price used a constant (mostly 0.5) for all the variables, a
randomized where , =0.8; =0.0001 was used in Reference
(Kitayama, Arakawa, and Yamazaki 2011). Constraints are handled using a penalty formulation:

(4.33)

Algorithm
The algorithm used in (Kitayama, Arakawa, and Yamazaki 2011) is as follows:

• 0: Set , and . Initialize the iteration counter =1.

• 1: Generate the particles (parameter vectors) at random.

• 2: Apply the following procedure to all particles:

– 2-1: Particle , denoted by , selects three particles ( , , ) at random, where .

– 2-2: Mutation: A new particle is generated by the mutation: .

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 125
Optimization methods

– 2-3: Crossover: A new particle is generated by the crossover between and .

– 2-4: Selection: The objective function is evaluated at and and the particle is updated
according to the following criteria:

• 3:

• 4: If , go to 2, else terminate.

In addition to 11 analytical benchmark problems with varying between 2 and 10, Kitayama et al.
analyzed two structural constrained minimization examples (the optimal design of a spring ( = 3)
and a truss topology optimization ( = 28)). As in (Storn and Price 1997), DE/rand/1/bin was used.
The performance of the DE was compared to several other algorithms – Generalized Random Tun-
neling Algorithm (GRTA) (Kitayama and Yamazaki 2005), Particle Swarm Optimization (PSO) (Eberhart
and Kennedy 1995) (Schutte 2005), GA (Golberg 1989), and Distributed GA and Simulated Annealing)
– and was concluded to be competitive. although GRTA outperformed DE on the topology bench-
mark.

References
Aluffi-Pentini, F, V Parisi, and F Zirilli. 1985. “Global Optimization and Stochastic Differential Equations.”
Journal of Optimization Theory and Applications 47: 1–16.

Eberhart, Russell, and James Kennedy. 1995. “Particle Swarm Optimization.” In Proceedings of the
IEEE International Conference on Neural Networks, 4:1942–48. Citeseer.

Golberg, David E. 1989. “Genetic Algorithms in Search, Optimization, and Machine Learning.” Addion
Wesley 1989 (102): 36.

Ingber, Lester et al. 1993. “Adaptive Simulated Annealing (ASA).” Global Optimization C-Code, Caltech
Alumni Association, Pasadena, CA.

Kitayama, Satoshi, Masao Arakawa, and Koetsu Yamazaki. 2011. “Differential Evolution as the Global
Optimization Technique and Its Application to Structural Optimization.” Applied Soft Computing 11
(4): 3792–3803.

Kitayama, Satoshi, and Koetsu Yamazaki. 2005. “Generalized Random Tunneling Algorithm for
Continuous Design Variables.”

Mühlenbein, Heinz, and Dirk Schlierkamp-Voosen. 1993. “Predictive Models for the Breeder Genetic
Algorithm i. Continuous Parameter Optimization.” Evolutionary Computation 1 (1): 25–49.

Press, William H, Saul A Teukolsky, William T Vetterling, and Brian P Flannery. 1992. Numerical Recipes
in c. Cambridge university press Cambridge.

Schutte, Jaco Francois. 2005. Applications of Parallel Global Optimization to Mechanics Problems.
University of Florida.

Storn, Rainer, and Kenneth Price. 1997. “Differential Evolution-a Simple and Efficient Heuristic for
Global Optimization over Continuous Spaces.” Journal of Global Optimization 11 (4): 341.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
126 of ANSYS, Inc. and its subsidiaries and affiliates.
 Single-objective optimization

Voigt, Hans-Michael. 2005. “Soft Genetic Operators in Evolutionary Algorithms.” In Evolution and
Biocomputation: Computational Models of Evolution, 123–41. Springer.

4.2.6. Hybrid methods

[Link]. One-Click Optimization (OCO)

Available in: ModelCenter, optiSLang

Keywords— Hybrid optimization, Success factor, MOP assisted optimization

The One-Click Optimization (OCO) method is a hybrid and dynamic optimization approach. It effi-
ciently combines direct (high-fidelity model) and MOP (low fidelity model) assisted search strategies.
OCO is a general purpose optimizer that automatically and iteratively selects the most suitable
optimization methods exposed by optiSLang.

OCO considers the allowable optimization methods available in optiSLang such as NLPQL (Sec-
tion [Link] (p. 87)), simplex (Section [Link] (p. 92)), EA (Section [Link] (p. 110)), and so on. It tackles
the shortcomings of most static optimization heuristics by using a dynamic selection scheme that
starts by exploring the design space to study the statistical features of the response(s). An in-house
selection criterion called the "success factor" is introduced to compare the optimization algorithms
scheme. This success factor considers the number of evaluations, and the estimated improvement
in terms of objectives and feasibility, in order to rank optimization algorithms regarding there ex-
pected performance for the application at hand.

One of the most innovative aspect of OCO is its dynamic and adaptive nature. In fact, it allows de-
pending on the behavior of the algorithms, either the successive use of different algorithms, or the
use of one algorithm. For example, if an algorithm is performing well OCO keeps using it. On the
contrary, if another algorithm looks more promising, OCO switches dynamically to the expected
more promising algorithm. OCO can run several optimizers in parallel to combine global and local
search.

This dynamic switch system can help avoid local minima, arbitrary changes, and the overhead that
comes with a manual selection.

The most relevant algorithmic steps of the OCO are depicted in Figure 4.25 (p. 128).

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 127
Optimization methods

Figure 4.25: OCO workflow

4.3. Multi-objective optimization

Most real-world optimization problems include more than one objective and a large number of design
variables. This leads to the general formulation of the multi-objective optimization problem

(4.34)

with as the vector of design variables and the constraints of inequality and
equality . All solutions which comply with both constraints and variable bounds form the feasible
-dimensional design space. Each solution is assigned to a vector describing
one point of the -dimensional objective space. This is the major difference to the single-objective
optimization problem where there is only one objective. The mapping between design space and ob-
jective space is illustrated in Figure 4.26 (p. 129).

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
128 of ANSYS, Inc. and its subsidiaries and affiliates.
 Multi-objective optimization

Figure 4.26: Design space (left) and objective space (right)

Several procedures have been developed to solve a multi-objective optimization problem (Branke et
al. 2008). They are classified in non-interactive and interactive methods. In non-interactive methods the
multi-objective optimization problem is transferred into a singe-objective problem. This can be done
e. g. by choosing a preferred objective function and introducing the remaining objectives as constraints
(see Section 4.3.3 (p. 133))

(4.35)

Another possibility is to combine all objectives in a single function by using individual weights (see
Section 4.3.2 (p. 132)):

(4.36)

Both methods require good knowledge about the optimization potential with respect to each objective
function and a clear idea about the importance of the different objectives with respect to each other.

In early stages of the design process such information are often rare and a transformation of the multi-
objective task into a single objective task is not possible a priori. In such cases interactive methods are
more promising. In optiSLang Pareto optimization as one of the most popular interactive methods is
available.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 129
Optimization methods

Figure 4.27: Recommended flowchart for multi-objective optimization

In Figure 4.27 (p. 130) the recommended flow for a multi-objective optimization procedure is shown:
initially the inputs, possible objectives and constraint functions are defined. Sensitivity analysis is used
to detect unimportant parameters and to check the objective functions with respect to possible conflicts.
Afterwards, a multi-objective optimization is performed to determine the optimization potential within
the conflicting objectives and to derive suitable weighting factors for a following single-objective op-
timization. Finally this single-objective optimization determines an optimal design.

References
Branke, J., K. Deb, K. Miettinen, and R. Slowinski. 2008. Multiobjective Optimization. Springer.

4.3.1. Pareto optimization

[Link]. Pareto optimality

Solving the multi-objective optimization problem requires a criterion which regards the distinct
objectives. One possible criterion is the Pareto dominance which is formulated for two decision
vectors and in the design space as follows

• Solution dominates solution if it is feasible and better or equal in all objectives and
better in at least one objective.

• Solution is indifferent to solution if none of both solutions dominates the respective

other.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
130 of ANSYS, Inc. and its subsidiaries and affiliates.
 Multi-objective optimization

Figure 4.28: Pareto dominating (filled circles) and dominated designs (unfilled circles) including
the resulting Pareto frontier of two conflicting objectives. ( and dominate but are
indifferent to each other)

Using the dominance criterion allows for creating a partial order between different decision vectors
as illustrated in Figure 4.28 (p. 131). A solution is Pareto optimal if there is no decision vector that
would improve one objective without causing a deterioration in at least one other objective. In
other words a solution is Pareto optimal if it is not dominated by any other solution. The term
Pareto is named after Vilfredo Pareto, an Italian economist who used the concept in his studies of
economic efficiency and income distribution.

If all objectives are equally important and no preferences are made a priori, the dominance of a
solution is the only way to determine if it is better than others. As a result the non-dominated
subset out of the feasible set of solutions constitutes the Pareto set. The corresponding points in
the objective space are called Pareto frontier. The following requirements concerning the multi-
objective optimization can be formulated:

• Find a set of solutions close to the Pareto optimal solutions (convergence).

• Find solutions which are diverse enough to represent the whole Pareto frontier (diversity).

Although the optimization produces a set of Pareto optimal solutions, the user is interested in
getting only one solution. Additional preferences have to be introduced, which needs comprehensive
knowledge about the problem, to select a solution that meets the requirements.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 131
Optimization methods

[Link]. Analysis of conflicting objectives

In order to obtain different trade-off solutions there must be conflicting objectives. The results of
a preceding sensitivity analysis should be used to check the objective functions concerning possible
conflicts.

Figure 4.29: Conflicting objectives of the damped oscillator (maximum amplitude vs.
eigen-frequency) including cluster analysis of the Pareto optimal designs in the objective
space (left) and in the design space (right)

In case of the oscillator example the damped eigen-frequency is used as second objective function
which should be minimized. The anthill plots of the samples used for the sensitivity analysis and
design exploration are shown in Figure 4.29 (p. 132). The figure indicates that the minimization of
the maximum amplitude is in conflict with the minimization of the damped eigen-frequency. In
the anthill plot all Pareto dominant designs can be identified. Further information can be gained
by subdividing the Pareto optimal designs in the objective space in clusters which are investigated
with respect to their region in the design space. The figure indicates e.g. that the designs of cluster
2 are highly concentrated in the objective space but widely distributed in the design space. Further-
more, all Pareto optimal designs are close to the design space boundaries.

4.3.2. Weighted sum

The weighted sum method is one of the most widely used approaches to solving multi-objective
optimization problems. It combines all objectives in a single function by using individual weights .

Examples of weighted sum function:

(4.37)

where with

(4.38)

where with and the target of The weighted sum method is simple and easy
to use. For convex problems it guarantees to find solutions on the entire Pareto optimal set.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
132 of ANSYS, Inc. and its subsidiaries and affiliates.
 Multi-objective optimization

Limitations
• For mixed problems (min-max), user needs to convert all the objectives into one type.

• Uniformly distributed set of weights does not guarantee a uniformly distributed set of Pareto-op-
timal solutions.

• Two different set of weights not necessarily lead to two different Pareto- optimal solutions.

• Cannot find certain Pareto-optimal solutions in the case of a nonconvex objective space.

4.3.3. ε-Constraint
The -constraint method is a simple way to transfer a multi-objective optimization problem as a
single-objective optimization problem with additional constraint. User keeps one of the objectives,
and treats the remaining objective as constraints. The initial optimization problem (Equation
4.35 (p. 129)) becomes:

(4.39)

Different Pareto-optimal solutions can be found using different values. This method is applicable
to either convex or non-convex problems.

Limitations: The values have to be chosen carefully to lie between the minimum and maximum
values for each objective function.

4.3.4. Response Surface based methods

[Link]. Adaptive Metamodel of Optimal Prognosis (AMOP)

Available in: optiSLang

Keywords— Metamodel-based, adaptive, MOP

In the multi-objective AMOP approach, the user can select between two different adaption strategies.
The default strategy will run an internal Pareto optimization using an Evolutionary algorithm as
discussed in section Section [Link] (p. 136). Based on the optimal designs found in the Pareto
frontier, an update of the support points is selected using a space-filling criterion in the objective
space. This approach will lead to a wide spread of the Pareto optimal designs as shown in Figure
4.30 (p. 134).

The second approach calculates the product of the probability of a possible improvement of the
individual objective functions and considers constraint conditions in the same way as in the single-
objective case:

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 133
Optimization methods

(4.40)

The second approach is recommended for multi-objective optimization tasks with more than 3
objectives functions, where a best compromise between the objectives is easier to interpret than
a multi-dimensional Pareto frontier.

Figure 4.30: Adaptive MOP: Pareto frontier update for the damped oscillator example using
the default space-filling update criterion (left) and by using the best compromise criterion
(right)

Limitations: Due to the internal use of Kriging and MLS approximation, the AMOP algorithm for
multi-objective optimization is not advised for more than 10 input parameters.

[Link]. Adaptive Multiple-Objective

Available in: DesignXplorer, optiSLang

Keywords— Meta-model based

Adaptive Multiple-Objective (AMO) combines a Kriging response surface (see Section 3.1.4 (p. 26))
and Multi-Objective Genetic Algorithm (MOGA: see Section [Link] (p. 141)). It allows you to either
generate a new sample set or use an existing set, providing a more refined approach than the
Screening method. Except when necessary, the optimizer does not evaluate all design points. The
general optimization approach is the same as MOGA, but a Kriging response surface is used. Part
of the population is "simulated" by evaluations of the Kriging response surface. The Kriging error
predictor reduces the number of evaluations used in finding the first Pareto front solutions. AMO
supports multiple objectives and multiple constraints.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
134 of ANSYS, Inc. and its subsidiaries and affiliates.
 Multi-objective optimization

Figure 4.31: AMO workflow

The AMO algorithm proceeds in the following steps:

1. Initial population of MOGA: the initial population of MOGA is used to construct the Kriging re-
sponse surfaces.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 135
Optimization methods

2. Kriging generation: a Kriging response surface is created for each output, based on the first
population and then enhanced during simulation with the addition of new design points.

3. MOGA: MOGA is run to generate the next population.

4. Evaluate the population: the population is evaluated on the Kriging response surface.

5. Error check: Kriging error predictor is checked for each design of the new population. If the error
for a given design point is acceptable, the approximated value provided by Kriging is included
in the next population to be run through MOGA (return to Step 3). If the error is not acceptable,
the design point is promoted as new design point to evaluate an the real solver. The new design
points are used to improve the Kriging response surface (return to Step 2) and are included in
the next population to be run through MOGA (return to Step 3).

6. Convergence check: the optimization is validated for convergence. MOGA converges when the
maximum allowable Pareto percentage has been reached. When this happens, the process is
stopped. If the optimization is not converged, the process continues to the next step.

7. Stopping criteria: if the optimization has not converged, it is validated for fulfillment of the
stopping criteria. When the maximum number of iterations has been reached, the process is
stopped without having reached convergence. If the stopping criteria have not been met, the
MOGA algorithm is run again (return to Step 3).

8. Conclusion: Steps 2 through 7 are repeated in sequence until the optimization has converged
or the stopping criteria have been met. When either of these things occurs, the optimization
concludes.

Limitations: AMO does not support discrete parameters. It is limited to continuous parameters,
including those with manufacturable values (list of continuous values).

4.3.5. Natured inspired methods

[Link]. Evolutionary algorithms

Available in: ModelCenter, optiSLang

The application of evolutionary algorithms for solving multi-objective optimization problems has
been established over the past two decades. The parallel search for a set of Pareto optimal solutions
is the major advantage of this method. The optiSLang implementation is based on the Strength
Pareto Evolutionary Algorithm 2 (SPEA2) (Zitzler, Laumanns, and Thiele 2001) and is characterized
by the following principles:

• Elitism is applied by using an archive of non-dominated individuals.

• Fitness assignment is based on the dominance criterion.

• Preservation of population diversity is realized by density estimation.

The algorithm starts with the random or manual initialization of the first population of size followed
by the evaluation of all individuals. For the first generation (number of parents) individuals are
selected out of the population for reproduction. After the evaluation of all newborns both archive
individuals and offspring are assigned new fitness values. The best solutions form the archive of

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
136 of ANSYS, Inc. and its subsidiaries and affiliates.
 Multi-objective optimization

the next generation. The optimization stops if the maximum number of generations has been
reached or the archive stagnated. In contrast to the single-objective nature-inspired optimization
methods, fitness assignment of the multi-objective EA always regards the whole population, because
the criterion is based on the comparison of individuals. The fitness assignment method is a domin-
ance-based ranking which takes into account by how many individuals a solution is dominated and
how many individuals a solution dominates. To preserve the diversity of the Pareto frontier an ad-
ditional criterion is introduced to the fitness assignment procedure. The distance of an individual
to its -th nearest neighbor serves as an estimator of density. The intention is to prefer solutions
in less crowded regions of the objective space such as boundary solutions.

Figure 4.32: Pareto frontier for the damped oscillator (left) and Pareto optimal designs plotted
in the design space (right) obtained by evolutionary algorithms

In Figure 4.32 (p. 137) the resulting Pareto frontier of the multi-objective evolutionary algorithm is
shown for the damped oscillator example. The Pareto frontier shows a good agreement with the
Pareto optimal designs from Figure 4.29 (p. 132). The Pareto frontier can now be used to judge
about the importance of the different optimization goals and to choose a suitable optimal design.
The efficiency of the multi-objective optimization methods can be dramatically improved especially
for high-dimensional optimization problems by choosing a suitable start population instead of a
pure random generation. Selected designs of a sensitivity analysis or previous optimization runs
shall be used for this purpose.

References
Zitzler, E., M. Laumanns, and L. Thiele. 2001. "SPEA2: Improving the Strength Pareto Evolutionary
Algorithm." 103. Computer Engineering; Networks Laboratory (TIK), Swiss Federal Institute of Tech-
nology (ETH) Zurich.

[Link]. Darwin algorithm

Available in: ModelCenter

Keywords— Natures-inspired, Evolutionary Algorithm

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 137
Optimization methods

Darwin algorithm (see Section [Link] (p. 113)) can also be applied to multi-objective optimization
problems. When more than one objectives are defined, Darwin will search for a series of best designs
(Pareto set) using non-dominated design search.

[Link]. Particle Swarm Optimization

Available in: ModelCenter, optiSLang

The Particle Swarm Optimization algorithm can also be applied to multi-objective optimization
problems. As mentioned for single-objective PSO (see Section [Link] (p. 115)) optiSLang selects the
fittest design as global leader. For the next iteration step the swarm will move in this direction. In
multi-objective optimization the fitness evaluation is different from the single-objective approach.

Fitness assignment of the multi-objective PSO is similar to EA (Section [Link] (p. 136)).

[Link]. Non-dominated Sorting Genetic Algorithm II

Available in: LS-OPT, ModelCenter

This algorithm was developed by Prof. Kalyanmoy Deb and his students in 2000 (Deb et al. 2002).
This algorithm first tries to converge to the Pareto optimal front and then it spreads solutions to
get diversity on the Pareto optimal front. Since this algorithm uses a finite population size, there
may be a problem of Pareto drift. To avoid that problem, Goel et al.(Goel et al. 2007) proposed
maintaining an external archive.

The implementation of this archived NSGA-II is shown in Figure 4.33 (p. 139), and described as follows:

1. Randomly initialize the parent population (size ). Initialize an empty archive.

2. Evaluate the population i.e., compute constraints and objectives for each individual.

3. Rank the population using non-domination criteria. Also compute the crowding distance (this
distance finds the relative closeness of a solution to other solutions in the function space and
is used to differentiate between the solutions on same rank).

4. Employ genetic operators - selection, crossover & mutation - to create a child population.

5. Evaluate the child population.

6. Combine the parent and child populations, rank them, and compute the crowding distance.

7. Apply elitism (defined in a following section): Select best individuals from the combined
population. These individuals constitute the parent population in the next generation.

8. Add all rank = 1 solutions to the archive.

9. Update the archive by removing all dominated and duplicate solutions.

10. If the termination criterion is not met, go to step 4. Otherwise, report the candidate Pareto op-
timal set in the archive.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
138 of ANSYS, Inc. and its subsidiaries and affiliates.
 Multi-objective optimization

Figure 4.33: Elitist non-dominated sorting genetic algorithm (NSGA-II). The shaded blocks are
not the part of original NSGA-II but additions to avoid Pareto drift.

Elitism is applied to preserve the best individuals. The mechanism used by NSGA-II algorithm for
elitism is illustrated in Figure 4.34 (p. 140). After combining the child and parent populations, there
are 2 individuals. This combined pool of members is ranked using non-domination criterion
such that there are individuals with rank . The crowding distance of individuals with the same
rank is computed. Steps in selecting individuals are as follows:

1. Set , and number of empty slots .

2. If ,

• Copy all individuals with rank " " to the new parent population.

• Reduce the number of empty slots by : .

• Increment " ": .

• Return to Step 2.

3. If ,0.

• Sort the individuals with rank " " in decreasing order of crowding distance.

• Select individuals.

• Stop.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 139
Optimization methods

Figure 4.34: Elitism in NSGA-II.

To preserve diversity on the Pareto optimal front, NSGA-II uses a crowding distance operator. The
individuals with same rank are sorted in ascending order of function values. The crowding distance
is the sum of distances between immediate neighbors, such that in Figure 4.35 (p. 140), the crowding
distance of selected individual is " ". The individuals with only one neighbor are assigned a very
high crowding distance. Note: It is important to scale all functions such that they are of the same
order of magnitude otherwise the diversity preserving mechanism would not work properly.

Figure 4.35: Illustration of non-domination criterion, Pareto optimal set, and Pareto optimal
front.

NSGA II improved the original NSGA by adding following capabilities:

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
140 of ANSYS, Inc. and its subsidiaries and affiliates.
 Multi-objective optimization

• Faster sorting algorithm to compute Pareto set.

• Elitist approach is used to preserve best designs from previous generations.

• In multi-objective optimization, finding diverse designs is important. If a design is close to another

design, a process called "sharing" is used to penalize such points. NSGA II uses more robust
sharing technique that does not require control parameters.

References
Deb, Kalyanmoy, Amrit Pratap, Sameer Agarwal, and TAMT Meyarivan. 2002. "A Fast and Elitist
Multiobjective Genetic Algorithm: NSGA-II." IEEE Transactions on Evolutionary Computation 6 (2):
182–97.

Goel, Tushar, Rajkumar Vaidyanathan, Raphael T Haftka, Wei Shyy, Nestor V Queipo, and Kevin
Tucker. 2007. "Response Surface Approximation of Pareto Optimal Front in Multi-Objective Optim-
ization." Computer Methods in Applied Mechanics and Engineering 196 (4-6): 879–93.

[Link]. Multi-Objective Genetic Algorithms

Available in: DesignXplorer

The Multi-Objective Genetic Algorithm (MOGA) is a hybrid variant of the popular NSGA-II (Non-
dominated Sorted Genetic Algorithm-II) based on controlled elitism concepts. It supports all types
of input parameters. The Pareto ranking scheme is done by a fast, non-dominated sorting method
that is an order of magnitude faster than traditional Pareto ranking methods. The constraint handling
uses the same non-dominance principle as the objectives. Therefore, penalty functions and Lagrange
multipliers are not needed. This also ensures that the feasible solutions are always ranked higher
than the infeasible solutions.

The first Pareto front solutions are archived in a separate sample set internally and are distinct from
the evolving sample set. This ensures minimal disruption of Pareto front patterns already available
from earlier iterations. The basic steps of the MOGA algorithm are as follows:

1. Initialize the population.

2. Evaluate the initial population members (calculate the values of the objective function(s) and
constraints for each population member).

3. Loop until converged, or stopping criteria reached.

a. Perform crossover.

b. Perform mutation.

c. Evaluate the new population.

d. Assess the fitness of each member of the population so as to select members to be replaced
in next step.

e. Replace population with members selected to continue in the next generation.

f. Apply niche pressure to the population to generate a more even and uniform sampling by
encouraging differentiation along the Pareto frontier.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 141
Optimization methods

g. Test for convergence.

4.3.6. Hybrid methods

[Link]. One-Click Optimization (OCO)

Available in: ModelCenter, optiSLang

Keywords— Hybrid optimization, Success factor

The One-Click Optimization (OCO) method is a hybrid and dynamic optimization approach. It effi-
ciently combines direct (high-fidelity model) and MOP (low fidelity model) assisted search strategies.
OCO is a general purpose optimizer that automatically and iteratively selects the most suitable
optimization methods exposed by optiSLang.

The main difference between the OCO approach for multi-objective optimization and the OCO ap-
proach for single-objective optimization (see Section [Link] (p. 127)) is the computation of the success
factor. In multi-objective applications, the hypervolume indicator (see Figure 4.36 (p. 143)) is an in-
tegral part of the success factor, used to evaluate the success of an optimization algorithm and
rank optimization algorithms with a scalar metric.

4.3.7. Performance metrics

[Link]. Number of nondominated points

This is the number of solutions of all nondominated solutions at any generation. Usually a higher
number of nondominated points are achieved when convergence is good.

[Link]. Spread
The spread of the front is calculated as the diagonal of the largest hypercube in the function space
that encompassed all points. A large spread is desired to find diverse trade-off solutions. The spread
measure is derived using the extreme solutions making it susceptible to the presence of a few
isolated points that could artificially improve the spread metric. An equivalent criterion might be
the volume of such a hypercube.

[Link]. Standard deviation of crowding distance

This complimentary criterion (to the spread metric) detects the presence of poorly distributed
solutions by estimating how uniformly the points are distributed in the Pareto optimal set. This
metric is defined as,

(4.41)

where is the crowding distance of the solution in the function or variable space. The boundary
points are assigned a crowding distance of twice the distance to the nearest neighbor. A small
value of the uniformity measure is desired to achieve a good distribution of points.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
142 of ANSYS, Inc. and its subsidiaries and affiliates.
 Multi-objective optimization

[Link]. Min/Max of objectives

This represents the range of individual objectives. A wide range represents more choices for the
designer.

[Link]. Hypervolume
The hypervolume indicator was first proposed as a method for assessing multiobjective optimization
algorithms (Zitzler and Thiele 1999). A dominated hypervolume metric tries to simultaneously es-
timate the convergence and spread characteristics by computing the union of the volume between
the optimal solutions and a reference point. For practical purposes, the nadir point of all solutions
is used as the reference point. While all the above metrics are obtained on a single set of solutions,
the following performance metrics are obtained by comparing multiple sets of solutions. These
metrics are helpful in determining the convergence.

Figure 4.36: Hypervolume metric of two conflicting objectives ( and ) with as reference
point.

References
Zitzler, Eckart, and Lothar Thiele. 1999. "Multiobjective Evolutionary Algorithms: A Comparative
Case Study and the Strength Pareto Approach." IEEE Transactions on Evolutionary Computation 3 (4):
257–71.

[Link]. Number of common points

This is the number of solutions that exist in both sets and . A large number of common
points is indicative of the high quality of solutions. The set of common solutions is represented as,
(4.42)

and is the size of set . This is a particularly good metrics when a large generation interval
is used.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 143
Optimization methods

[Link]. Number of new nondominated solutions

This metrics denotes the number of nondominated solutions that were evolved during the current
generation interval. The set of such solutions is represented as,
(4.43)

A large number of new solutions relative to the total archive size indicates that the new solutions
are still being evolved and hence convergence is not yet achieved.

[Link]. Number of old dominated solutions n(Q)

This metrics denotes the number of nondominated solutions in the older archive that were
dominated by the solutions in the current archive . The set of dominated solutions is,
(4.44)

A large number of new solutions relative to the total archive size indicates that the new solutions
are still being evolved and hence convergence is not yet achieved.

[Link]. Consolidation ratio

This represents the fraction of archive that has evolved up to the generation. This is com-
puted as the ratio of the number of members in archive that are also present in the archive
(non-dominated solutions) to the size of archive . Mathematically,

(4.45)

This metric represents the proportion of potentially converged solutions in the archive. In the early
phase of a multi-objective evolutionary algorithm (MOEA) simulation, a large fraction of the non-
dominated solutions in the archive would be dominated by the solutions in archive due to
evolution, thus resulting in a small fraction of surviving solutions i.e., small value of the consolidation
ratio. However, significantly better solutions evolve in the later phases such that a large proportion
of solutions in the archive remain non-dominated with respect to new solutions leading to a
high consolidation ratio. In the limiting case, the consolidation ratio approaches one.

[Link]. Improvement ratio

This represents the fraction of archive dominated by the new solutions in archive . This is
computed as the ratio of the number of members in archive that are dominated by the solutions
in archive (dominated solutions) to the size of archive . Mathematically,

(4.46)

The archive includes all non-dominated members of archive so no member of the archive
is dominated. The improvement ratio quantifies the extent of improvement in the quality of
evolved solutions. This metric has a high value in the early phase of simulation which gradually
converges to zero when convergence is achieved.

More information about these performance metrics can be obtained from (Goel and Stander 2010).

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
144 of ANSYS, Inc. and its subsidiaries and affiliates.
 Optimization Glossary

References
Goel, Tushar, and Nielen Stander. 2010. "A Non-Dominance-Based Online Stopping Criterion for
Multi-Objective Evolutionary Algorithms." International Journal for Numerical Methods in Engineering
84 (6): 661–84.

4.4. Optimization Glossary

Adaptive Response Surface (ARSM). Response surface based method to solve single-objective optim-
ization. The objective and constraint functions are approximated by polynomial regression. This approach
generates a new adapted design scheme in each iteration step.

ARSM. See Adaptive Response Surface Method.

Adaptive Multi-Objective. Response surface based method to solve multi-objective optimization. Ge-
netic algorithm coupling solver evaluations and response surface evaluations.

Adaptive Simulated Annealing. Global stochastic optimization algorithm that mimics the metallurgical
annealing process.

Adaptive Single-Objective. Response surface based method to solve single-objective optimization. It

employs advanced refinement methods to provide the global optima.

AMO. See Adaptive Multi-Objective.

ASA. See Adaptive Simulated Annealing.

ASO. See Adaptive Single-Objective.

Constraint. An absolute limit on a response variable specified in terms of an upper or lower limit.

Constrained optimization. The mathematical optimization of a function subject to specified limits on

other functions.

Covariance Matrix Adaptation. Stochastic, or randomized, evolution strategy for real-parameter

(continuous domain) optimization of non-linear, non-convex functions.

CMA. See Covariance Matrix Adaptation.

Design of Experiments. See experimental design.

Design parameter. See design variable.

Design formula. A simple mathematical expression which gives the response of a design when the
design variables are substituted. See response surface.

Design space. A region in the -dimensional space of the design variables ( through ) to which
the design is limited. The design space is specified by upper and lower bounds on the design variables.
Response variables can also be used to bound the design space.

Design surface. The response variable as a function of the design variables, used to construct the for-
mulation of a design problem. (See also response surface, design rule).

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 145
Optimization methods

Design sensitivity. The gradient vector of the response. The derivatives of the response function in
terms of the design variables. df /dxi.

Design variable. An independent design parameter which is allowed to vary in order to change the
design. Symbolized by ( or (vector containing several design variables)).

Differential Evolution. Evolutionary Algorithm (EA) where its mutation operator is quite different, using
a geometric approach that is motivated by the moves performed in the Downhill simplex method.

Domain reduction. The reduction of the region of interest in the design space during the optimization
process.

Downhill simplex. Also known as the Nelder-Mead method. Direct search method to find the minimum
(or maximum) of a multivariate function.

D-optimal. The state of an experimental design in which the determinant of the moment matrix of the
least squares formulation is maximized.

DOE. Design of Experiments. See experimental design.

DSA. Design sensitivity analysis.

Efficient Global Optimization. Bayesian global optimization (BO) technique that consists in sampling
the point that maximizes the so-called expected improvement (EI).

Evolutionary Algorithm (EA). Stochastic search methods that mimic processes of natural biological
evolution.

Experimental Design. The selection of designs to enable the construction of a design response surface.
Sometimes referred to as the Point Selection Scheme.

EA. See Evolutionary Algorithm.

EGO. See Efficient Global Optimization.

Feasible Design. A design which complies with the constraint bounds.

Function. A mathematical expression for a response variable in terms of design variables. Often used
interchangeably with "response". Symbolized by f.

Functionally efficient. See Pareto optimal.

Function evaluation. Using a solver to analyze a single design and produce a result. See Simulation.

Global approximation. A design function which is representative of the entire design space.

Global Optimization. The mathematical procedure for finding the global optimum in the design space.
e.g. Genetic Algorithm, Particle Swarm, etc.

Global Sensitivity Analysis. A sensitivity analysis method which uses Sobol indices.

Gradient vector. A vector consisting of the derivatives of a function f in terms of a number of variables
x1 to xn. s = [df /dxi]. See Design Sensitivity.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
146 of ANSYS, Inc. and its subsidiaries and affiliates.
 Optimization Glossary

GSA. See Global Sensitivity Analysis. History. Response history containing two columns of (usually time)
data generated by a simulation.

Importance. See Weight.

Infeasible Design. A design which does not comply with the constraint functions. An entire design
space or region of interest can sometimes be infeasible.

Iteration. A cycle involving an experimental design, function evaluations of the designs, approximation
and optimization of the approximate problem.

Kriging. A Metamodeling technique using Bayesian regression.

Latin Hypercube Sampling. The use of a constrained random experimental design as a point selection
scheme for response approximation.

Leapfrog Optimization. Gradient-based method to solves constrained optimization problems.

Local Approximation. See Gradient vector.

LFOPC. Leapfrog optimizer for constrained minimization. See Leapfrog Optimization.

MDO. Multidisciplinary design optimization.

MOEA. Multi-objective evolutionary algorithm.

Metamodeling. The construction of surrogate design models such as polynomial response surfaces,
Artificial Neural Networks or Kriging surfaces from simulations at a set of design points.

Multidisciplinary design optimization (MDO). The inclusion of multiple disciplines in the design op-
timization process. In general, only some design variables need to be shared between the disciplines
to provide limited coupling in the optimization of a multidisciplinary target or objective.

Multi-objective. An objective function which is constituted of more than one objective. Symbolized
by F.

Multi-objective Optimization (MOO). Multi-objective optimization is the procedure for constructing

a Pareto optimal front.

Multi-criteria. Refers to optimization problems in which several criteria are considered.

MISQP. Solves mixed-integer nonlinear programming problems by a modified sequential quadratic

programming (SQP) method. See Sequential quadratic programming.

MOO. Multi-objective Optimization.

MOP. Metamodel of Optimal Prognosis. General framework for an automatic model testing and selection
in optiSLang.

MP. Mathematical Programming. Mathematical optimization.

Noise. See random error.

NLPQLP. Solves non-linear programming problems by a sequential quadratic programming (SQP)

method. See Sequential quadratic programming.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 147
Optimization methods

Objective. A function of the design variables that the designer wishes to minimize or maximize. If there
exists more than one objective, the objectives have to be combined mathematically into a single ob-
jective. Symbolized by .

Optimal design. The methodology of using mathematical optimization tools to improve a design iter-
atively with the objective of finding the "best" design in terms of predetermined criteria.

Optimization strategy. A strategy for metamodel-based optimization such as Single Stage, Sequential
or Sequential with Domain Reduction.

Pareto Front. Set of designs that are all Pareto-optimal. See Pareto optimal.

Pareto optimal. A multi-objective design is Pareto-optimal if none of the objectives can be improved
without at least one objective being affected adversely. A Pareto optimal front can be constructed using
optimization.

Particle Swarm Optimization. Nature-inspired optimization method that imitates the social behavior
of a swarm.

Preference function. A function of objectives used to combine several objectives into a single one
suitable for the standard MP formulation. Preprocessor.

Process simulation. The use of computer programming, computer vision, and feedback to simulate
manufacturing techniques.

PI-BO. Probabilistic Inference for Bayesian Optimization. It is an algorithm based on the so-called
Bayesian Optimization with DIM-GP (Deep Infinite Mixture Gaussian Process) as internal response surface.

PSO. See Particle Swarm Optimization.

Region of interest. A sub-region of the design space. Usually defined by a mid-point design and a
range of each design variable. Usually dynamic.

Response. A numerical indicator of the performance of the design. A function of the design variables
approximated using a metamodel which can be used for optimization. Symbolized by f. Collected over
all design iterations for plotting.

Response quantity. See response.

Response Surface. A mathematical expression which relates the response variables to the design
parameters. Typically computed using statistical methods.

Response variable. A response which is substituted as an input variable in a downstream stage of a

multi-stage process. Response-expressions are allowed.

RSM. Response Surface Methodology.

Sampling. In general, Sampling is synonymous with Point Selection or Experimental Design.

Sensitivity. See Design sensitivity.

Sequential quadratic programming. SQP is a class of algorithms for solving non-linear optimization
problems (NLP). SQP methods solve a sequence of optimization subproblems, each of which optimizes
a quadratic model of the objective subject to a linearization of the constraints.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
148 of ANSYS, Inc. and its subsidiaries and affiliates.
 Optimization Glossary

Simulation. The analysis of a physical process or entity in order to compute useful responses. See
Function evaluation.

Solver. A computational tool used to analyze a structure or fluid using a mathematical model. The ex-
ecutable software for such a tool. See Discipline.

Stage. A distinct step or operation in a process which typically reads input, processes the input and
produces a result. Example: run a solver. Different stages can be dependent on one another.

Stochastic. Involving or containing random variables. Involving probability or chance.

Stochastic Design Improvement. local, single-objective optimization procedure that improves a pro-
posed design by using a simple stochastic approach. In each iteration step, the best design is taken as
new center point for the sampling of the next. The ranges of the sampling scheme are adapted in each
iteration step.

Stopping Criterion. A mathematical criterion for terminating an iterative procedure.

SDI. See Stochastic Design Improvement.

SQP. See Sequential quadratic programming.

Target. A desired value for a response. The optimizer will not use this value as a rigid constraint. Instead,
it will try to get as close as possible to the specified value.

Tolerance optimization. Optimization which includes tolerances on the input variables.

Trade-off curve. A curve constructed using Pareto optimal designs.

Weight. A measure of importance of a response function or objective. Typically varies between 0 and
1.

Weighted Sum Method. Transfers a multi-objective optimization problem into a singe-objective

problem. It combines all objectives in a single function by using individual weights.

Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
of ANSYS, Inc. and its subsidiaries and affiliates. 149
 Release 2025 R2 - © ANSYS, Inc. All rights reserved. - Contains proprietary and confidential information
150 of ANSYS, Inc. and its subsidiaries and affiliates.