Appendix: Test Problems

D.1 Unconstrained Problems¶

D.1.1 Slanted Quadratic Function¶

This is a smooth two-dimensional function suitable for a first test of a gradient-based optimizer:

where $\beta \in [0, 2)$. A $\beta$ value of zero corresponds to perfectly circular contours. As $\beta$ increases, the contours become increasingly slanted. For $\beta=2$, the quadratic becomes semidefinite, and there is a line of weak minima. For $\beta > 2$, the quadratic is indefinite, and there is no minimum. An intermediate value of $\beta=3/2$ is suitable for first tests and yields the contours shown in Figure D.1.
Global minimum: $f(x^*) = 0$ at $x^* = (0, 0)$.

D.1.2 Rosenbrock Function¶

The two-dimensional Rosenbrock function, shown in Figure D.2, is also known as Rosenbrock’s valley or banana function. This function was introduced by Rosenbrock,1 who used it as a benchmark problem for optimization algorithms.

The function is defined as follows:

This became a classic benchmarking function because of its narrow turning valley. The large difference between the maximum and minimum curvatures, and the fact that the principal curvature directions change along the valley, makes it a good test for quasi-Newton methods.

The Rosenbrock function can be extended to $n$ dimensions by defining the sum,

Global minimum: $f(x^*) = 0.0$ at $x^* = (1,1,\ldots,1)$.

Local minimum: For $n \ge 4$, a local minimum exists near $x=(-1,1,\ldots,1)$.

D.1.3 Bean Function¶

The “bean” function was developed in this book as a milder version of the Rosenbrock function: it has the same curved valley as the Rosenbrock function but without the extreme variations in curvature. The function, shown in Figure D.3, is

Global minimum: $f(x^*) = 0.09194$ at $x^* = (1.21314, 0.82414)$.

D.1.4 Jones Function¶

This is a fourth-order smooth multimodal function that is useful to test global search algorithms and also gradient-based algorithms starting from different points. There are saddle points, maxima, and minima, with one global minimum. This function, shown in Figure D.4 along with the local and global minima, is

Global minimum: $f(x^*) = -13.5320$ at $x^* = (2.6732, -0.6759)$.
Local minima: $f(x) = -9.7770$ at $x = (-0.4495, 2.2928)$.
$f(x) = -9.0312$ at $x = (2.4239, 1.9219)$.

D.1.5 Hartmann Function¶

The Hartmann function is a three-dimensional smooth function with multiple local minima:

where

A slice of the function, at the optimal value of $x_1 = 0.1148$, is shown in Figure D.5.
Global minimum: $f(x^*) = -3.86278$ at $x^* = (0.11480, 0.55566, 0.85254)$.

D.1.6 Aircraft Wing Design¶

We want to optimize the rectangular planform wing of a general aviation-sized aircraft by changing its wingspan and chord (see Example 1.1). In general, we would add many more design variables to a problem like this, but we are limiting it to a simple two-dimensional problem so that we can easily visualize the results.

The objective is to minimize the required power, thereby taking into account drag and propulsive efficiency, which are speed dependent. The following describes a basic performance estimation methodology for a low-speed aircraft. Implementing it may not seem like it has much to do with optimization. The physics is not important for our purposes, but practice translating equations and concepts into code is an important element of formulating optimization problems in general.

In level flight, the aircraft must generate enough lift to equal the required weight, so

We assume that the total weight consists of a fixed aircraft and payload weight $W_0$ and a component of the weight that depends on the wing area $S$—that is,

The wing can produce a certain lift coefficient ($C_L$) and so we must make the wing area ($S$) large enough to produce sufficient lift. Using the definition of lift coefficient, the total lift can be computed as

where $q$ is the dynamic pressure and

If we use a rectangular wing, then the wing area can be computed from the wingspan ($b$) and the chord ($c$) as

The aircraft drag consists of two components: viscous drag and induced drag. The viscous drag can be approximated as

For a fully turbulent boundary layer, the skin friction coefficient, $C_f$, can be approximated as

In this equation, the Reynolds number is based on the wing chord and is defined as follows:

where $\rho$ is the air density, and $\mu$ is the air dynamic viscosity.

The form factor, $k$, accounts for the effects of pressure drag. The wetted area, $S_\text{wet}$, is the area over which the skin friction drag acts, which is a little more than twice the planform area. We will use

The induced drag is defined as

where $e$ is the Oswald efficiency factor. The total drag is the sum of induced and viscous drag, $D = D_i + D_f$.

Our objective function, the power required by the motor for level flight, is

where $\eta$ is the propulsive efficiency. We assume that our electric propellers have a Gaussian efficiency curve (real efficiency curves are not Gaussian, but this is simple and will be sufficient for our purposes):

In this problem, the lift coefficient is provided. Therefore, to satisfy the lift requirement in (8), we need to compute the velocity using (11) and (10) as

This is the same problem that was presented in Example 1.2 of Chapter 1. The optimal wingspan and chord are $b = 25.48$ m and $c = 0.50$ m, respectively, given the parameters. The contour and the optimal wing shape are shown in Figure D.6.

Because there are no structural considerations in this problem, the resulting wing has a higher wing aspect ratio than is realistic. This emphasizes the importance of carefully selecting the objective and including all relevant constraints.

The parameters for this problem are given as follows:

D.1.7 Brachistochrone Problem¶

The brachistochrone problem is a classic problem proposed by Johann Bernoulli (see Section 2.2 for the historical background). Although this was originally solved analytically, we discretize the model and solve the problem using numerical optimization. This is a useful problem for benchmarking because you can change the number of dimensions.

A bead is set on a wire that defines a path that we can shape. The bead starts at some $y$-position $h$ with zero velocity. For convenience, we define the starting point at $x=0$.

From the law of conservation of energy, we can then find the velocity of the bead at any other location. The initial potential energy is converted to kinetic energy, potential energy, and dissipative work from friction acting along the path length, yielding the following:

Now that we know the speed of the bead as a function of $x$, we can compute the time it takes to traverse an differential element of length $\mathrm{d} s$:

To discretize this problem, we can divide the path into linear segments. As an example, Figure D.7 shows the wire divided into four linear segments (five nodes) as an approximation of a continuous wire. The slope $s_i = (\Delta y/\Delta x)_i$ is then a constant along a given segment, and $y(x) = y_i + s_i (x - x_i)$. Making these substitutions results in

Performing the integration and simplifying (many steps omitted here) results in

where $\Delta x_i = (x_{i+1} - x_i)$, and $\Delta y_i = (y_{i+1} - y_i)$. The objective of the optimization is to minimize the total travel time, so we need to sum up the travel time across all of our linear segments:

Minimization is unaffected by multiplying by a constant, so we can remove the multiplicative constant for simplicity (we see that the magnitude of the acceleration of gravity has no effect on the optimal path):

The design variables are the $n-2$ positions of the path parameterized by $y_i$. The endpoints must be fixed; otherwise, the problem is ill-defined, which is why there are $n-2$ design variables instead of $n$. Note that $x$ is a parameter, meaning that it is fixed. We could space the $x_i$ any reasonable way and still find the same underlying optimal curve, but it is easiest to just use uniform spacing. As the dimensionality of the problem increases, the solution becomes more challenging. We will use the following specifications:

• Starting point: $(x, y) = (0, 1)$ m.

• Ending point: $(x, y) = (1, 0)$ m.

• Kinetic coefficient of friction $\mu_k = 0.3.$

The analytic solution for the case with friction is more difficult to derive, but the analytic solution for the frictionless case ($\mu_k = 0$) with our starting and ending points is as follows:

where $a = 0.572917$ and $\theta \in [0, 2.412]$.

D.1.8 Spring System¶

Consider a connected spring system of two springs with lengths of $l_1$ and $l_2$ and stiffnesses of $k_1$ and $k_2$, fixed at the walls as shown in Figure D.8. An object with mass $m$ is suspended between the two springs. It will naturally deform such that the sum of the gravitational and spring potential energy, $E_p$, is at the minimum.

The total potential energy of the spring system is

where $\Delta l_1$ and $\Delta l_2$ are the changes in length for the two springs. With respect to the original lengths, and displacements $x_1$ and $x_2$ as shown, they are defined as

This can be minimized to determine the final location of the object.

With initial lengths of $l_1 = 12$ cm, $l_2 = 8$ cm; spring stiffnesses of $k_1 = 1.0$ N$\cdot$cm, $k_2 = 10.0$ N$\cdot$cm; and a force due to gravity of $mg = 7\text{N}$, the minimum potential energy is at $(x_1, x_2) = (2.7852, 6.8996)$. The contour of $E_k$ with respect to $x_1$ and $x_2$ is shown in Figure D.9.

The analytic derivatives can also be computed for use in a gradient-based optimization. The derivative of $E_p$ with respect to $x_1$ is

where the partial derivatives of $\Delta l_1$ and $\Delta l_2$ are

By letting $\mathcal{L}_1 = \sqrt{\left(l_1 + x_1\right)^2 + x_2^2}$ and $\mathcal{L}_2 = \sqrt{\left(l_2 - x_1\right)^2 + x_2^2}$, the partial derivative of $E_p$ with respect to $x_1$ can be written as

Similarly, the partial derivative of $E_p$ with respect to $x_2$ can be written as

D.2 Constrained Problems¶

D.2.1 Barnes Problem¶

The Barnes problem was devised in a master’s thesis2 and has been used in various optimization demonstration studies. It is a good starter problem because it only has two dimensions for easy visualization while also including constraints.

The objective function contains the following coefficients:

For convenience, we define the following quantities:

The objective function is then:

There are three constraints of the form $g(x) \le 0$:

The problem also has bound constraints. The original formulation is bounded from $[0, 80]$ in both dimensions, in which case the global optimum occurs in the corner at $x^* = [80, 80]$, with a local minimum in the middle. However, for our usage, we preferred the global optimum not to be in the corner and so set the bounds to $[0, 65]$ in both dimensions. The contour of this function is plotted in Figure D.10.
Global minimum: $f(x^*) = -31.6368$ at $x^* = (49.5263, 19.6228)$.
Local minimum: $f(x) = -17.7754$ at $x = (65, 65)$.

D.2.2 Ten-Bar Truss¶

The 10-bar truss is a classic optimization problem.3 In this problem, we want to find the optimal cross-sectional areas for the 10-bar truss shown in Figure D.11. A simple truss finite-element code set up for this particular configuration is available in the book code repository.

The function takes in an array of cross-sectional areas and returns the total mass and an array of stresses for each truss member.

The objective of the optimization is to minimize the mass of the structure, subject to the constraints that every segment does not yield in compression or tension. The yield stress of all elements is $25 \times 10^3$ psi, except for member 9, which uses a stronger alloy with a yield stress of $75 \times 10^3$ psi. Mathematically, the constraint is

where the absolute value is needed to handle tension and compression (with the same yield strength for tension and compression). Absolute values are not differentiable at zero and should be avoided in gradient-based optimization if possible. Thus, we should put this in a mathematically equivalent form that avoids absolute value.

Each element should have a cross-sectional area of at least 0.1 in$^2$ for manufacturing reasons (bound constraint). When solving this optimization problem, you may need to scale the objective and constraints.

Although not needed to solve the problem, an overview of the equations is provided. A truss element is the simplest type of structural finite element and only has an axial degree of freedom. The theory and derivation for truss elements are simple, but for our purposes, we skip to the result. Given a two-dimensional element oriented arbitrarily in space (Figure D.12), we can relate the displacements at the nodes to the forces at the nodes through a stiffness relationship.

In matrix form, the equation for a given element is $K_{e} d = q$. In detail, the equation is

where the displacement vector is $d = [u_1, v_1, u_2, v_2]$. The meanings of the variables in the equation are described in Table 1.

The stress in the truss element can be computed from the equation $\sigma = S_e d$, where $\sigma$ is a scalar, $d$ is the same vector as before, and the element $S_e$ matrix (really a row vector because stress is one-dimensional for truss elements) is

The global structure (an assembly of multiple finite elements) has the same equations, $K d = q$ and $\sigma = S d$, but now $d$ contains displacements for all of the nodes in the structure, $d = [x_1, x_2, \ldots, x_n]$. If we have $n$ nodes and $m$ elements, then $q$ and $d$ are $2n$-vectors, $K$ is a $(2n \times 2n)$ matrix, $S$ is an $(m \times 2n)$ matrix, and $\sigma$ is an $m$-vector. The elemental stiffness and stress matrices are first computed and then assembled into the global matrices. This is straightforward because the displacements and forces of the individual elements add linearly.

After we assemble the global matrices, we must remove any degrees of freedom that are structurally rigid (already known to have zero displacement). Otherwise, the problem is ill-defined, and the stiffness matrix will be ill-conditioned.

Given the geometry, materials, and external loading, we can populate the stiffness matrix and force vector. We can then solve for the unknown displacements from

With the solved displacements, we can compute the stress in each element using