optimization.tex

% Copyright (C) 2013  Jim Turner
%
% This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License.
% To view a copy of this license, visit https://creativecommons.org/licenses/by-sa/4.0/.

\documentclass{article}
\title{Engineering Design Optimization}

\input{preamble.tex}
\setcounter{tocdepth}{2}

\begin{document}

\input{firstpage.tex}

\section{Standard Form}\label{sec:standard-form}

Find \(\vec{x} =
\begin{Bmatrix}
  x_1 \\
  x_2 \\
  \vdots \\
  x_n
\end{Bmatrix}\)
which minimizes \(f(\vec{x})\) subject to the constraints
\(\begin{aligned}
  g_j(\vec{x}) &\le 0, \ j = 1, 2, \dots, m \\
  h_k(\vec{x}) &\le 0, \ k = 1, 2, \dots, l
\end{aligned}\)

\section{Conditions for Optimality}\label{sec:conditions-optimality}

\subsection{Unconstrained}\label{sec:unconstrained-conditions}

\begin{description*}
  \item[Necessary conditions]~
    \begin{itemize*}
      \item \(\del F(x) = 0\)
      \item \(H(x)\) is positive semi-definite or positive definite
    \end{itemize*}
  \item[Sufficient conditions]~
    \begin{itemize*}
      \item \(H(x)\) is positive definite
    \end{itemize*}
\end{description*}

\subsection{Constrained}\label{sec:constrained-conditions}

Kuhn--Tucker Conditions are necessary for a relative minimum at \(\vec{x}^*\) and are sufficient to
ensure a global minimum at \(\vec{x}^*\) for convex programming problems:
\begin{enumerate*}
  \item \(\vec{x}^*\) is feasible (meets constraints)
  \item \(\lambda_j g_j(\vec{x}^*) = 0\) where \(j = 1, 2, \dots, m\) and \(\lambda_j \ge 0\)
  \item \(0 = \del f(\vec{x}^*) + \sum_{j=1}^m{\lambda_j \del g_j(\vec{x}^*)} +
    \sum_{k=1}^l{\lambda_{m+k} \del h_k(\vec{x}^*)}\)
    where \(\lambda_{m+k}\) are unrestricted in sign
\end{enumerate*}
To solve for \(\lambda_j\) for problems with only inequality constraints, let
\[\vec{B} = \del f(\vec{x}^*) \textrm{ and } A =
\begin{bmatrix}
  \del g_1(\vec{x}^*) & \del g_2(\vec{x}^*) & \cdots & \del g_n(\vec{x}^*)
\end{bmatrix}\]
where \(A\) only contains active constraints, then
\[\vec{\lambda} = -{\left[\transpose{A} A\right]}^{-1} \transpose{A} \vec{B}\]

\section{1-D Optimization}\label{sec:1-d-optimization}

\subsection{Bounding\slash{}Range Finding}\label{sec:bounding}

\begin{enumerate*}
\item \label{itm:choose-init-and-step} Choose \(x^0\) and a step size \(\Delta\). Set \(k = 0\). The
  larger the value of \(\Delta\), the fewer the function calls required for bounding but the larger
  the final bounds will be.
\item If \(f(x^0 - \abs{\Delta}) \ge f(x^0) \ge f(x^0 + \abs{\Delta})\),
  then \(\Delta\) is positive \\
  else if \(f(x^0 - \abs{\Delta}) \le f(x^0) \le f(x^0 + \abs{\Delta})\),
  then \(\Delta\) is negative \\
  else go to step~\ref{itm:choose-init-and-step}.
\item \label{itm:calc-next-x} \(x^{k+1} = x^k + 2^k \Delta\)
\item If \(f(x^{k+1}) < f(x^k)\),
  then increment \(k\) and go to step~\ref{itm:calc-next-x} \\
  else the minimum is between \((x^{k-1}, x^{k+1})\)
\end{enumerate*}

\subsection{Interpolation\slash{}Approximation Methods}\label{sec:interp-approx}

Interpolation\slash{}approximation methods approximate the objective function with simpler functions
for which the minima are known, then iteratively improve their approximations. They are more
unpredictable and less robust than region elimination methods but converge faster for continuous and
well-conditioned problems.

\subsubsection{Three-Point Quadratic Non-Derivative Optimization with Refinement}\label{sec:quadratic-interpolation}

\begin{enumerate*}
\item Start with three points: \((x_1, f_1)\), \((x_2, f_2)\), \((x_3, f_3)\) where \(x_1\),
  \(x_3\) bound the minimum and \(x_2 \in (x_1, x_3)\). Typically, \(x_2 = \frac12(x_1 + x_3)\).
\item \label{itm:define-interp} Define an interpolation function \(\tilde{f}(x) = a_0 + a_1(x - x_1)
  + a_2(x - x_1)(x - x_2) \approx f(x)\)
\item Find \(f_\mrm{min} = \min{(f_1, f_2, f_3)}\) and associated \(x_\mrm{min}\).
\item To find the interpolation function, \(\tilde{f}(x)\):
\begin{align*}
  a_0 &= f_1 \\
  a_1 &= \frac{f_2 - f_1}{x_2 - x_1} \\
  a_2 &= \frac{1}{x_3 - x_2}\left(\frac{f_3 - f_1}{x_3 - x_1} - \frac{f_2 - f_1}{x_2 - x_1}\right)
\end{align*}
Then at the minimum of \(\tilde{f}(x)\), \(\left.\pd{\tilde{f}}{x}\right|_{\tilde{x}^*} = a_1 +
a_2(\tilde{x}^* - x_2) + a_2(\tilde{x}^* - x_1) = 0\). Solving for \(\tilde{x}^*\):
\begin{align*}
\tilde{x}^* &= \frac{x_2 + x_1}{2} - \frac{a_1}{2a_2} \\
\tilde{f}^* &= f(\tilde{x}^*)
\end{align*}
\item Check convergence and if converged, stop:
\[\abs{\frac{\tilde{x}^* - x_\mrm{min}}{x_\mrm{min}}} \le \epsilon_x \qquad \qquad
\abs{\frac{\tilde{f}^* - f_\mrm{min}}{f_\mrm{min}}} \le \epsilon_f\]
\item Save the best point (either \(\tilde{x}^*\) or \(x_\mrm{min}\)) and the two points that
  bracket it. Relabel the saved points as \(x_1\), \(x_2\), \(x_3\). Recalculate \(f_1\), \(f_2\),
  \(f_3\). Go to step~\ref{itm:define-interp}.
\end{enumerate*}

\subsubsection{Newton--Raphson Method}\label{sec:newton-raphson}

This method uses first and second derivative information to speed up convergence. However,
discontinuities are a problem, for some problems it can be expensive to get the derivative, and the
algorithm can diverge in some cases.

\begin{enumerate*}
\item Select a value \(x^0\).
\item \label{itm:newton-linear-approx} Build a linear approximation of \(f'\):
\[\tilde{f}'(x) = f'(x^k) + f''(x^k)(x - x^k)\]
\item Solve \(\tilde{f}'(x^{k+1}) = 0\) for \(x^{k+1}\):
\[x^{k+1} = x^k - \frac{f'(x^k)}{f''(x^k)}\]
\item Check for convergence (when \(\abs{x^{k+1} - x^k}\) is small). If not converged, go to
  step~\ref{itm:newton-linear-approx}.
\end{enumerate*}

\subsubsection{Two-Point Cubic Optimization with Refinement}\label{sec:cubic-interpolation}

This interpolation-based method uses first derivative information to generate splines for
\(\tilde{f}(x)\) and \(\tilde{f}'(x)\).

\subsection{Region Elimination Methods}\label{sec:region-elimination}

Region elimination methods iteratively eliminate subintervals of the design space from
consideration. They generally eliminate the same percentage of the space on each iteration. They are
more robust than interpolation\slash{}approximation method, particularly for discontinuous or
ill-conditioned functions, but may be slower for well-conditioned problems.

\subsubsection{Interval Halving Method}\label{sec:interval-halving}

After \(n\) calls to \(f\), the space is reduced to about \({\left(\frac12\right)}^{n/2}\) of its
original size.

\begin{enumerate*}
\setcounter{enumi}{-1}
\item Should already know \((x_\mrm{L}, f_\mrm{L})\) and \((x_\mrm{R}, f_\mrm{R})\) from bounding
  algorithm.
\item Calculate the following:
  \begin{align}
    x_\mrm{m} &= \frac{x_\mrm{L} + x_\mrm{R}}{2} \nonumber \\
    f_\mrm{m} &= f(x_\mrm{m}) \nonumber \\
    L &= x_\mrm{R} - x_\mrm{L} \label{eq:interval-halving-L}
  \end{align}
\item \label{itm:update-x1-x2} Calculate the following:
  \begin{align*}
    x_1 &= x_\mrm{L} + \frac{L}{4} \\
    f_1 &= f(x_1) \\
    x_2 &= x_\mrm{R} - \frac{L}{4} \\
    f_2 &= f(x_2)
  \end{align*}
\item Set \(x_\mrm{L}\), \(x_\mrm{m}\), and \(x_\mrm{R}\) to the three points that make a V shape.
\item Compute \(L\) using equation~\ref{eq:interval-halving-L} and check convergence. Go to
  step~\ref{itm:update-x1-x2} if not converged.
\end{enumerate*}

\subsubsection{Golden Section Method}\label{sec:golden-section-region}

Define the golden ratio conjugate as \(\Phi = \frac{\sqrt5-1}{2} \approx 0.61803\). After \(n\)
calls to \(f\), the space is reduced to about \(\Phi^n\) of its original size.

\begin{enumerate*}
\setcounter{enumi}{-1}
\item Should already know \((x_\mrm{L}, f_\mrm{L})\) and \((x_\mrm{R}, f_\mrm{R})\) from bounding
  algorithm.
\item Calculate the following:
  \begin{align}
    x_1 &= \Phi x_\mrm{L} + (1-\Phi) x_\mrm{R} \label{eq:golden-section-x1} \\
    f_1 &= f(x_1) \label{eq:golden-section-f1} \\
    x_2 &= (1-\Phi) x_\mrm{L} + \Phi x_\mrm{R} \label{eq:golden-section-x2} \\
    f_2 &= f(x_2) \label{eq:golden-section-f2}
  \end{align}
\item \label{itm:update-golden-section} If \(f_2 > f_1\), then
  \begin{itemize*}
  \item \(x_\mrm{R} := x_2\), \(f_\mrm{R} := f_2\)
  \item \(x_2 := x_1\), \(f_2 := f_1\)
  \item Compute the new values of \(x_1\), \(f_1\) using
    equations~\ref{eq:golden-section-x1}~and~\ref{eq:golden-section-f1}.
  \end{itemize*}
  else
  \begin{itemize*}
  \item \(x_\mrm{L} := x_1\), \(f_\mrm{L} := f_1\)
  \item \(x_1 := x_2\), \(f_1 := f_2\)
  \item Compute the new values of \(x_2\), \(f_2\) using
    equations~\ref{eq:golden-section-x2}~and~\ref{eq:golden-section-f2}.
  \end{itemize*}
\item Compute \(\epsilon = \frac{x_\mrm{R} - x_\mrm{L}}{L_0}\) and check convergence. If not
  converged, go to step~\ref{itm:update-golden-section}.
\end{enumerate*}

\subsubsection{Bisection Method}\label{sec:bisection-region}

This method uses first derivative information to eliminate half of the space on each iteration.

\begin{enumerate*}
\setcounter{enumi}{-1}
\item Should already know \((x_\mrm{L}, f_\mrm{L})\) and \((x_\mrm{R}, f_\mrm{R})\) from bounding
  algorithm. Also compute \(f'_\mrm{L} < 0\) and \(f'_\mrm{R} > 0\).
\item \label{itm:update-bisection-xm} Compute \(x_\mrm{m} = \frac{x_\mrm{L} + x_\mrm{R}}{2}\).
\item Compute \(f_\mrm{m} = f(x_\mrm{m})\) and \(f'_\mrm{m} = f'(x_\mrm{m})\).
\item If \(f'_\mrm{m} > 0\), then \(x_\mrm{R} := x_\mrm{m}\), else \(x_\mrm{L} := x_\mrm{m}\).
\item Compute \(\epsilon = \frac{x_\mrm{R} - x_\mrm{L}}{L_0}\) and check convergence. If not
  converged, go to step~\ref{itm:update-bisection-xm}.
\end{enumerate*}

\subsection{Hybrid Methods}

Often, the best general method is to use a region elimination method for a few iterations to reduce
the size of the bounds, then use an interpolation\slash{}approximation method to quickly converge on
the minimum.

\section{\textit{n}-D Optimization}

Unconstrained \(n\)-D optimization methods are typically formulated as a series of 1-D searches. The
individual methods determine the direction \(\vec{S}^q\) for each search. The update relation is
typically written as the following for iteration \(q\):
\[\vec{x}^{q+1} = \vec{x}^q + \alpha^* \vec{S}^q\]
where the 1-D optimization varies \(\alpha^*\) to find the minimum of \(F(\vec{x}^q + \alpha^*
\vec{S}^q)\).

\subsection{Zero Order Methods}

\subsubsection{Brute Force}

Don't actually use this method unless your objective function is very fast to evaluate.

\begin{enumerate*}
\item Pick a base point \(\vec{x}^0\) and set \(q = 0\).
\item \label{itm:brute-force-eval} Evaluate \(F(\vec{x}^q)\).
\item Pick sample points around the base point and evaluate \(F\) at those points.
\item Set the new base point \(\vec{x}^{q+1}\) to the point with the lowest value of \(F\). Stop if
  converged or increment \(q\) and go to step~\ref{itm:brute-force-eval}.
\end{enumerate*}

\subsubsection{Global Random Search}

This can be useful to get an understanding of your search space. Otherwise, it's not a good idea in
general.

\begin{enumerate*}
\item \label{itm:random-pick} Pick a random \(\vec{x}^q\) and evaluate \(F(\vec{x}^q)\).
\item If converged (\(n\) or \(F\) cutoff), stop. Otherwise, go to step~\ref{itm:random-pick}.
\end{enumerate*}

\subsubsection{Powell's Method of Conjugate Directions}

\paragraph{Advantages}
\begin{itemize}
\item Faster than a series of univariate moves because it doesn't slow down as much near the
  optimum. (The conjugate search direction generally points toward the optimum.)
\item Minimizes a quadratic function in a finite number of steps.
\end{itemize}

\paragraph{Disadvantages}
\begin{itemize}
\item If a search gains no improvement, conjugacy is lost and the method breaks down.
\item The method slows down near the optimum. If this becomes too extreme, it may be helpful to
  restart the algorithm from the latest point.
\end{itemize}

\paragraph{Methods}
\begin{enumerate*}
\item Set the initial search directions to the identity matrix
  \(S^0
  = \begin{bmatrix}\vec{S}^0_1 & \vec{S}^0_2 & \cdots & \vec{S}^0_n\end{bmatrix}
  = I_n\) and set \(q = 0\).
\item \label{itm:powell-n-searches} For each \(\vec{S}^q_i\) in \(S^q\):
  \begin{enumerate*}
  \item Find the corresponding \(\alpha^*\) to minimize \(F(\vec{x}^q_i + \alpha^* \vec{S}^q_i)\).
  \item Set \(\vec{x}^q_{i+1} = \vec{x}^q_i + \alpha^* \vec{S}^q_i\).
  \end{enumerate*}
\item Calculate the conjugate direction \(\vec{S}^q_\mathrm{c} = \vec{x}^q_n - \vec{x}^q_0\).
\item Find the corresponding \(\alpha^*\) to minimize \(F(\vec{x}^q_n + \alpha^*
  \vec{S}^q_\mathrm{c})\).
\item Set \(\vec{x}^{q+1}_0 = \vec{x}^q_n + \alpha^* \vec{S}^q_\mathrm{c}\).
\item Set the new search matrix \(S^{q+1} =
\begin{bmatrix}
\vec{S}^q_2 & \vec{S}^q_3 & \cdots & \vec{S}^q_n & \vec{S}^q_\mathrm{c}
\end{bmatrix}\).
\item If not converged, increment \(q\) and go to step~\ref{itm:powell-n-searches}.
\end{enumerate*}

\subsection{First Order Methods}

First order methods are usually more efficient than zero order methods, but they need gradient
information and perform poorly when the gradient is not continuous.

\subsubsection{Steepest Descent (Cauchy's Method)}

This method is simple to implement but slows down near the optimum. The method is to always set the
search direction to \(\vec{S}^q = -\del F(\vec{x}^q)\).

\subsubsection{Fletcher--Reeves Conjugate Direction Method}

This method uses steepest descent as a first move, and then uses a conjugate version of successive
gradients for subsequent moves.

\begin{enumerate*}
\item Pick a starting point \(\vec{x}^0\) and set \(q = 0\).
\item Set \(\vec{S}^0 = -\del F(\vec{x}^0)\).
\item Find \(\alpha^*\) and set \(\vec{x}^1 = \vec{x}^0 + \alpha^* \vec{S}^0\). Increment \(q\).
\item \label{itm:fletcher-beta} Set \(\displaystyle \beta^{q-1} = \frac{{\abs{\del F(\vec{x}^q)}}^2}{{\abs{\del F(\vec{x}^{q-1})}}^2}\).
\item Set \(\vec{S}^q = -\del F(\vec{x}^q) + \beta^{q-1} \vec{S}^{q-1}\).
\item Find \(\alpha^*\) and set \(\vec{x}^{q+1} = \vec{x}^q + \alpha^* \vec{S}^q\).
\item Check convergence. If not converged, go to step~\ref{itm:fletcher-beta}.
\end{enumerate*}

\subsection{Newton's Method}

The search direction for Newton's Method is:
\[\vec{S}^q = -{\left[H(\vec{x}^q)\right]}^{-1} \del F(\vec{x}^q)\]
The basic method is to set \(\alpha^* = 1\). As long as \(H(\vec{x}^q)\) is positive definite, this
is a good move. The method can diverge or overshoot the minimum. To reduce overshoots, a modified
update relation is to solve for the \(\alpha^*\) that gives the 1-D minimum for \(F(\vec{x}^q +
\alpha^* \vec{S}^q)\) instead of just setting it to \(1\).

\subsection{Variable Metric Methods (Quasi-Newton)}

These methods iteratively improve their approximation of the Hessian or its inverse. This allows
them to work similarly to Newton's method, but with only first order information.

\subsubsection{Davidon--Fletcher--Powell (DFP) Method}

This method approximates the inverse of the Hessian.

\begin{enumerate*}
\item Set
  \begin{itemize*}
  \item the initial point \(\vec{x}^0\)
  \item the inverse Hessian approximation \(H^0 = I_n\)
  \item the vector \(\vec{c}^0 = \del F(\vec{x}^0)\)
  \item the convergence criterion \(\varepsilon\)
  \item the iteration number \(q = 0\)
  \end{itemize*}
\item \label{itm:dfp-converge} If \(\norm{\vec{c}^q} < \varepsilon\), stop because the method has converged. Otherwise,
  continue to the next step.
\item Set \(\vec{S}^q = -H^q \vec{c}^q\).
\item Find \(\alpha^*\) to minimize \(F(\vec{x}^q + \alpha^* \vec{S}^q)\).
\item Set \(\vec{x}^{q+1} = \vec{x}^q + \alpha^* \vec{S}^q\).
\item Set \(H^{q+1} = H^q + B^q + C^q\) where:
  \begin{align*}
    B^q &= \frac{\vec{p}\transpose{\vec{p}}}{\transpose{\vec{p}}\vec{y}} \\
    C^q &= \frac{-\vec{z}\transpose{\vec{z}}}{\transpose{\vec{y}}\vec{z}} \\
    \vec{p} &= \vec{x}^{q+1} - \vec{x}^q \\
    \vec{y} &= \del F(\vec{x}^{q+1}) - \del F(\vec{x}^q) \\
    \vec{z} &= H^q \vec{y}
  \end{align*}
\item Increment \(q\) and go to step~\ref{itm:dfp-converge}.
\end{enumerate*}

\subsubsection{Broyden--Fletcher--Goldfarb--Shanno (BFGS) Method}

This method approximates the Hessian directly.

\begin{enumerate*}
\item Set
  \begin{itemize*}
  \item the initial point \(\vec{x}^0\)
  \item the Hessian approximation \(H^0 = I_n\)
  \item the vector \(\vec{c}^0 = \del F(\vec{x}^0)\)
  \item the convergence criterion \(\varepsilon\)
  \item the iteration number \(q = 0\)
  \end{itemize*}
\item \label{itm:bfgs-converge} If \(\norm{\vec{c}^q} < \varepsilon\), stop because the method has converged. Otherwise,
  continue to the next step.
\item Solve \(H^q\vec{S}^q = -\vec{c}^q\) for \(\vec{S}^q\).
\item Find \(\alpha^*\) to minimize \(F(\vec{x}^q + \alpha^* \vec{S}^q)\).
\item Set \(\vec{x}^{q+1} = \vec{x}^q + \alpha^* \vec{S}^q\).
\item Set \(H^{q+1} = H^q + D^q + E^q\) where:
  \begin{align*}
    D^q &= \frac{\vec{y}\transpose{\vec{y}}}{\transpose{\vec{y}}\vec{p}} \\
    E^q &= \frac{\vec{c}\transpose{\vec{c}}}{\transpose{\vec{c}}\vec{S}} \\
    \vec{p} &= \vec{x}^{q+1} - \vec{x}^q \\
    \vec{y} &= \del F(\vec{x}^{q+1}) - \del F(\vec{x}^q)
  \end{align*}
\item Increment \(q\) and go to step~\ref{itm:bfgs-converge}.
\end{enumerate*}

\section{Scaling\slash{}Normalization of Design Variables and Constraints}

\subsection{Scaling Design Variables}

Scaling design variables helps some of the algorithms perform more effectively. One method is to use
elements from the Hessian matrix:
\[\tilde{\vec{x}} = D \vec{x} \quad \textrm{where} \quad D =
\begin{bmatrix}
  \frac{1}{\sqrt{H_{11}}} & 0 & \cdots & 0 \\
  0 & \frac{1}{\sqrt{H_{22}}} & \cdots & 0 \\
  \vdots & \vdots & \ddots & \vdots \\
  0 & 0 & \cdots & \frac{1}{\sqrt{H_{nn}}}
\end{bmatrix}\]
Another method is to \(\vec{x}\) directly:
\[\tilde{\vec{x}} = D \vec{x} \quad \textrm{where} \quad D =
\begin{bmatrix}
  \frac{1}{\abs{x_1}} & 0 & \cdots & 0 \\
  0 & \frac{1}{\abs{x_2}} & \cdots & 0 \\
  \vdots & \vdots & \ddots & \vdots \\
  0 & 0 & \cdots & \frac{1}{\abs{x_n}}
\end{bmatrix}\]

\subsection{Normalization of Constraints}

Normalizing constraints provides similar benefits to normalizing design variables. To normalize,
they can be written as:
\[g_\textrm{normalized}(\vec{x}) = \frac{g(\vec{x})}{g_\textrm{allowable}} - 1 \le 0\]

\subsection{Eliminating Bounds}

Lower and upper bounds on \(x_i\) can often be eliminated by expressing \(x_i\) in terms of another
variable \(y_i\). For example, if \(\ell_i \le x_i \le u_i\), then two options are:
\begin{align*}
  x_i &= \ell_i + (u_i - \ell_i) \sin^2 y_i \\
  x_i &= \ell_i + (u_i - \ell_i) \frac{y_i^2}{1 + y_i^2}
\end{align*}

\section{Sequential Unconstrained Minimization Techniques (SUMT)}

These techniques transform a constrained optimization problem into an unconstrained one. They
typically use an unconstrained pseudo-objective function for the optimization. This is
pseudo-objective function is typically written as \(\Phi(\vec{x}, r_\mrm{p})\) where \(r_\mrm{p}\)
is a penalty parameter and \(P(\vec{x})\) is a penalty function:
\[\Phi(\vec{x}, r_\mrm{p}) = F(\vec{x}) + r_\mrm{p} P(\vec{x})\]

\subsection{Exterior Penalty Function}

This method approaches the optimum starting from the infeasible region.

\paragraph{Advantages}
\begin{itemize}
\item Finding an initial point is easy.
\end{itemize}

\paragraph{Disadvantages}
\begin{itemize}
\item The method never gets back into the feasible region (although it gets close).
\end{itemize}

\paragraph{Method}
The penalty function is defined as follows:
\[P(\vec{x})
= \sum_{j=1}^m{\left[\max\left(0, g_j(\vec{x})\right)\right]}^2
+ \sum_{k=1}^\ell{\left[h_k(\vec{x})\right]}^2\]
This function has a slope of 0 at the constraint boundary and has a continuous slope for
\(\Phi\). However, the second derivative is not continuous. The penalty parameter \(r_\mrm{p}\) is
chosen small initially and is increased as the method progresses, typically by the relation:
\[r_\mrm{p}^\textrm{next} = \gamma r_\mrm{p} \quad \textrm{where usually} \quad \gamma \in [3, 10]\]

\subsection{Interior Penalty Function}

This method approaches the optimum starting from the feasible region.

\paragraph{Advantages}
\begin{itemize}
\item If you stop the method early, you still have a feasible design.
\end{itemize}

\paragraph{Disadvantages}
\begin{itemize}
\item The pseudo-objective function is discontinuous at the boundaries of the constraints.
\item The constraints need to be normalized.
\item Region elimination (which can be slower) is safer for the 1-D search.
\end{itemize}

\paragraph{Method}
The pseudo-objective function is defined as follows:
\[\Phi(\vec{x}, r_\mrm{p}, r'_\mrm{p})
= F(\vec{x}) + r'_\mrm{p} \sum_{j=1}^m{\frac{-1}{g_j(\vec{x})}}
+ r_\mrm{p} \sum_{k=1}^\ell{\left[h_k(\vec{x})\right]}^2\]
\(r_\mrm{p}\) behaves the same as for the exterior penalty function. \(r'_\mrm{p}\) starts large and then decreases:
\begin{align*}
  {r'_\mrm{p}}^\textrm{init} &= (0.1 \textrm{ to } 1) \left(\frac{f(\vec{x})}{-\sum_{j=1}^m{\frac{1}{g_j(\vec{x})}}}\right) \\
  {r'_\mrm{p}}^\textrm{next} &= \gamma r'_\mrm{p} \quad \textrm{where} \quad \gamma \in \{0.1, 0.2, 0.5\}
\end{align*}

\paragraph{Extended Method}
This extension helps keep the penalty function from ``blowing up'' near the optimum. The portion
of \(P(\vec{x})\) for each \(g(\vec{x})\) is defined as:
\[\textrm{portion of } P(\vec{x}) \textrm{ for each } g(\vec{x}) =
\begin{cases}
  \frac{-1}{g_j(\vec{x})} & \textrm{if } g_j(\vec{x}) < \varepsilon \\
  \frac{-2\varepsilon - g_j(\vec{x})}{\varepsilon^2} & \textrm{if } g_j(\vec{x}) > \varepsilon
\end{cases}\] where \(\varepsilon\) is a small negative number defined as
\(\varepsilon = -c {(r'_\mrm{p})}^a\)

\subsection{Augmented Lagrange Multiplier (ALM) Method}

This method is usually superior to the interior and exterior penalty function
methods. This section needs to be completed.

\section{Simplex Method}

\begin{description*}
\item[Canonical form]
  Minimize \(F(\vec{x}) = c_1x_1 + c_2x_2 + c_3x_3 + \cdots + c_nx_n\) subject to
  \begin{align*}
    a_{11}x_1 + a_{12}x_2 + a_{13}x_3 + \cdots a_{1n}x_n &= b_1 \\
    a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + \cdots a_{2n}x_n &= b_2 \\
    \vdots \\
    a_{n1}x_1 + a_{n2}x_2 + a_{n3}x_3 + \cdots a_{nn}x_n &= b_n \\
    x_1,x_2,\cdots,x_n &\ge 0
  \end{align*}
  where all \(b \ge 0\).
\end{description*}

Things to do to turn nonlinear standard form into simplex standard form:
\begin{enumerate}
\item If an \(x_i\) is not constrained to be \(\ge 0\), then substitute \(x_i = x_i' - x_i''\) where \(x_i',x_i'' \ge 0\).
\item Use slack variables to turn inequality constraints into equalities:
  \begin{align*}
    a_1x_1 + a_2x_2 + \cdots \le b &\implies a_1x_1 + a_2x_2 + \cdots + s_1 = b \\
    a_1x_1 + a_2x_2 + \cdots \ge b &\implies a_1x_1 + a_2x_2 + \cdots - s_1 = b
  \end{align*}
\item Add artificial variables to constraints that were initially \(=\) or \(\ge\), e.g.:
  \[a_1x_1 + a_2x_2 + \cdots \ge b \implies a_1x_1 + a_2x_2 + \cdots - s_1 + a_1 = b\]
\end{enumerate}

The tableau is constructed from the simplex standard form. Variables with a 1 in
one coefficient and a 0 everywhere else are in the basis, and are assumed to be
nonzero. All other variables are nonbasic and are assumed to be zero. The steps
in the simplex method are:
\begin{enumerate}
\item Get artificial variables into the basis, then work to eliminate the \(w\)s, bottom-up. If a \(w\) cannot be eliminated, no feasible solution can be found.
\item Figure out which variable \(j\) you want to bring into the basis to reduce \(F\)
  the fastest. This is the variable with the most negative coefficient in the last row of the tableau.
\item Pivot about the row \(i\) that has the smallest value of \(b_i / a_{ij}\) for rows where \(a_{ij} > 0\).
\item The method has converged when all of the coefficients in the last row are
  \(\ge 0\) or the values are negative and grow unbounded. When converged:
  \begin{itemize*}
  \item If all of the coefficients of non-basic variables in the last row are \(\ge 0\), the solution is unique.
  \item If some of the coefficients of non-basic variables in the last row are \(= 0\), the solution is not unique.
  \item If some of the coefficients of non-basic variables in the last row are \(< 0\), the solution is unbounded.
  \end{itemize*}
\end{enumerate}

\section{Sequential Linear Programming}

\subsection*{Steps}

\begin{enumerate*}
\item Find initial point
\item \label{itm:linearize} Linearize about point
\item Apply move limits (20--25\% of DV range initially, about 5\% near optimum)
\item Minimize linear problem
\item Check convergence; if not converged, go to step~\ref{itm:linearize}.
\end{enumerate*}

\subsection*{Notes}

\begin{itemize*}
\item Advantages
  \begin{itemize*}
  \item practical
  \item easy to implement
  \item fast for lots of DVs
  \end{itemize*}
\item Disadvantages
  \begin{itemize*}
  \item have to select move limits and how to shrink them
  \end{itemize*}
\end{itemize*}

\section{Method of Feasible Directions}

This method is only designed to handle inequality constraints. It uses a pushoff factor to avoid violating constraints. A method is
\begin{description*}
\item[Usable if] \(\transpose{\del F(\vec{x})} \cdot \vec{S} \le 0\)
\item[Feasible if] \(\transpose{\del g_j(\vec{x})} \cdot \vec{S} \le 0\)
\end{description*}

TODO

\section{Generalized Reduced Gradient (GRG) Method}

Find \(\vec{S}\) such that any active constraint remains precisely active for some small move.

\begin{enumerate*}
\item Convert the problem into standard GRG form. There are \(n\)
  design variables, and all inequality constraints are transformed into equality
  constraints by adding \(m\) slack variables:

  Minimize \(F(\vec{x})\) such that
  \begin{align*}
    g_j(\vec{x}) + x_{j+n} &= 0 & j &= 1,m \\
    h_k(\vec{x}) &= 0 & k &= 1,\ell \\
    x_i^L \le x_i &\le x_i^U & i &= 1,n \\
    x_{j+n} &\ge 0 & j &= 1,m
  \end{align*}

  Each equality constraint makes one independent variable a dependent variable. The variables are split into two parts:
  \[\vec{x} = \begin{Bmatrix}\vec{y} \\ \vec{z}\end{Bmatrix} \quad
  \begin{array}{l} n-\ell \text{ independent variables} \\ m + l \text{ dependent variables} \end{array}\]

\item Find the starting point \(\vec{x}^0\).
\item \label{itm:reduced-grad} Calculate the search direction from the reduced gradient:
  \[\vec{S} = -\od{F(\vec{x})}{\vec{y}} = -\left[\del_{\vec{y}}F - \transpose{\left({[D]}^{-1} [C]\right)} \del_{\vec{z}} F\right]\]
  where
  \begin{align*}
    [C] &= \begin{bmatrix}
      \pd{h_1}{y_1} & \cdots & \pd{h_1}{y_{n-\ell}} \\
      \vdots & \ddots & \vdots \\
      \pd{h_{m+\ell}}{y_1} & \cdots & \pd{h_{m+\ell}}{y_{n - \ell}} \\
    \end{bmatrix} \\
    [D] &= \begin{bmatrix}
      \pd{h_1}{z_1} & \cdots & \pd{h_1}{z_{m+\ell}} \\
      \vdots & \ddots & \vdots \\
      \pd{h_{m+\ell}}{z_1} & \cdots & \pd{h_{m+\ell}}{z_{m+\ell}} \\
    \end{bmatrix}
  \end{align*}
\item Find the \(\alpha^*\) that minimizes \(F\) along the \(\vec{S}\) search direction, and calculate the corresponding \(\vec{y}\).
\item Find the \(\vec{z}\) corresponding to \(\vec{y}\) by iterating until \(\dif{\vec{z}} = 0\).
  \begin{align*}
    \dif{\vec{z}} &= -{[D]}^{-1}\left(\vec{h}(\vec{x}) + [C]\dif{\vec{y}}\right) \\
    \vec{z}^{k+1} &= \vec{z}^k + \dif{\vec{z}}
  \end{align*}
\item The \(\vec{y}\) and \(\vec{z}\) vectors give the new design vector \(\vec{x}\).
  If not converged, go to step~\ref{itm:reduced-grad}.
\end{enumerate*}

\section{Simulated Annealing}

This method simulates the thermal annealing of critically heated solids. At high
temperatures, the atoms move more freely. As the temperature is reduced, the
atoms move less freely and are more ordered, minimizing the internal energy. The
process depends on the cooling rate.

\begin{enumerate*}
\item Start with \(\vec{x}^0\) and a high value for \(T\). A good starting value
  for \(T\) is the average of some randomly selected \(F(\vec{x})\). Select a
  maximum number of iterations \(n\) for this temperature. A good value for \(n \in [50, 100]\).
\item \label{itm:random-design-proposal} Randomly generate a proposed design in the vicinity of the current point.
\item The probability of the proposed design \(\vec{x}^{i+1}\) depends on \(\Delta E\):
\begin{align*}
  \Delta E &= E^{i+1} - E^{i} = F(\vec{x}^{i+1}) - F(\vec{x}^i) \\
  P(E^{i+1}) &= \min\left(1, \mrm{e}^{-\Delta E/(kT)}\right)
\end{align*}
Note that if \(\Delta E < 0\), the proposed point is always accepted. If \(\Delta E > 0\),
the point is accepted less often when \(T\) is small. Increment \(i\).
\item Check for convergence. Can stop if \(\Delta f\) is small enough or \(T\) is small enough.
\item If \(i < n\), go to step~\ref{itm:random-design-proposal}. Otherwise, let
  \(T = cT\) and \(i = 1\), then go to step~\ref{itm:random-design-proposal}. A
  good value of \(c \in [0.4, 0.6]\). More complex cooling schedules are in the literature.
\end{enumerate*}

\subsection*{Features}

\begin{itemize*}
\item The quality of the solution is not influenced by the starting point,
  although the computational effort could increase.
\item It does not need to be continuous or differentiable.
\item Convergence is not influenced by convexity of the design space.
\item \(\vec{x}^i\) are not required to be positive.
\item Can solve mixed-integer, discrete, or continuous problems.
\item To account for constraints, use a penalty function.
\end{itemize*}

\section{Genetic Algorithms}

These methods are modeled after the natural process of evolution. They have the following steps:
\begin{enumerate*}
\item Select initial population
\item \label{itm:ga-selection} Select parents
\item Crossover
\item Recombine
\item Mutation and go to step~\ref{itm:ga-selection}
\end{enumerate*}

\subsection*{Design Encoding Methods}

\begin{itemize*}
\item Vector of integer
\item Vector of real numbers
\item Vector of binary bits (can represent integer vectors this way)
  \begin{itemize*}
  \item Flexible for many different data types
  \item Natural representation for computers
  \item Easy to crossover
  \item Can be unnatural to encode real numbers
  \end{itemize*}
\item Data structure
\end{itemize*}

\subsection*{Initial Population Generation}

\begin{itemize*}
\item Random
\item Pattern (e.g. latin hypercube)
\end{itemize*}

\subsection*{Evaluating Design Fitness}

\begin{itemize*}
\item Must be transitive (\(F(a) \ge F(b) \ge F(c) \implies F(a) \ge F(c)\))
\item The fitness must be increasingly better closer to the optimum
\item Can be equal to the objective function
\item Creative fitness functions can improve performance
\end{itemize*}

\subsection*{Selecting Parents and Producing New Designs}

\subsubsection*{Selection}

\begin{itemize*}
\item Roulette wheel (fitness proportionate selection)
  \begin{tabular}{ccc}
    \toprule
    Design & Fitness & Probability / \% \\
    \midrule
    1 & \(F_1\) & \(F_1 / \sum{F}\) \\
    2 & \(F_2\) & \(F_2 / \sum{F}\) \\
    3 & \(F_3\) & \(F_3 / \sum{F}\) \\
    \(\vdots\) & \(\vdots\) & \(\vdots\) \\
    \midrule
    Total & \(\sum{F}\) & 100\% \\
    \bottomrule
  \end{tabular}
\item Tournament
  \begin{itemize*}
  \item \(b\) strings are chosen at random
  \item most fit individual in group wins and becomes parent
  \item \(n\) tournaments to find \(n\) parents
  \end{itemize*}
\item Truncation: top \(1/q\) individuals get \(q\)
  copies in the mating pool, or select only the top \% of parents
\end{itemize*}

\subsubsection*{Crossover}

\begin{itemize*}
\item Strategies to generate new population
  \begin{itemize*}
  \item Use a probability for crossover, then children replace parents
  \item Select a set of parents to create a fixed number of children, which are
    added back to the population. Then, downsize the population.
  \end{itemize*}
\item Binary representations
  \begin{itemize*}
  \item \(k\)-point crossover strategy for binary strings
    \begin{itemize*}
    \item \(1 \le k \le \ell - 1\) is defined by the user, where \(\ell\) is the number of bits in a string
    \item randomly choose \(k\) splits in the string and swap every other section to generate two children
    \end{itemize*}
  \item uniform crossover strategy for binary strings: every bit can be swapped based on some probability (usually 0.5)
  \end{itemize*}
\item Continuous representations
  \begin{itemize*}
  \item Arithmetic crossover: the children are the following, where \(\lambda\) is generated randomly
    \begin{align*}
      \text{child 1} &= \lambda(\text{parent 1}) + (1 - \lambda)(\text{parent 2}) \\
      \text{child 2} &= (1 - \lambda)(\text{parent 1}) + \lambda(\text{parent 2})
    \end{align*}
  \item BLX-\(\alpha\) generates one child. For each DV,
    \begin{align*}
      P_\mrm{min} &= \min(\text{parent 1 DV},\ \text{parent 2 DV}) \\
      P_\mrm{max} &= \max(\text{parent 1 DV},\ \text{parent 2 DV}) \\
      \text{child DV} &\in \left[P_\mrm{min} - I\cdot\alpha, P_\mrm{max} + I\cdot\alpha\right]
    \end{align*}
  \item Linear crossover generates three children, and the best two of the three are selected.
    \begin{align*}
      C_1 &= \tfrac12 P_1 + \tfrac12 P_2 \\
      C_2 &= \tfrac32 P_1 - \tfrac12 P_2 \\
      C_3 &= -\tfrac12 P_1 + \tfrac32 P_2
    \end{align*}
  \end{itemize*}
\end{itemize*}

\subsubsection*{Mutation}

\begin{itemize*}
\item Binary: bitwise not
\item Real values
  \begin{itemize*}
  \item mutation over distribution (e.g. normal)
  \item random reassignment within bounds
  \end{itemize*}
\end{itemize*}

\subsection*{Implementation Issues}

\begin{itemize*}
\item handling constraints: could use penalty function
\item convergence
  \begin{itemize*}
  \item average fitness of population has stagnated
  \item homogeneous population
  \item best design is not changing
  \item maximum number of generations
  \item maximum number of \(F\) calls
  \end{itemize*}
\end{itemize*}

\subsection*{Establish Parameters for Algorithm}

\begin{itemize*}
\item Population size: generally 5 to 10 times the number of DVs (limit 200)
\item Crossover rate/type: typically mate 75\% to 80\% of pairs
\item Mutation rate/type: typically 0.5\% for bits and 3--5\% for real values
\end{itemize*}

\section{Particle Swarm Optimization}

Particles' state is their position and velocity. They remember the best location
they have discovered so far and communicate good positions to each other. On
each iteration, they adjust their positions and velocities.

\subsection*{Rules}

\begin{description*}
\item[Cohesion] Stick together
\item[Separation] Don't come too close to each other
\item[Alignment] Follow general heading of the flock
\end{description*}

\subsection*{Principles}

\begin{itemize*}
\item When a particle finds good information, it transmits this to others.
\item All other particles gravitate toward the target, but not directly.
\item Each particle has independent thinking and past memory.
\end{itemize*}

\subsection*{Method}

Assume unconstrained (with some bounds on DVs).
\begin{enumerate*}
\item Generate initial population (typically 20--30 particles). Typically, set the initial velocities to 0.
\item \label{itm:swarm-update} For \(i\)th iteration, find the following for each particle \(j\):
  \[\vec{v}_j(i) = \vec{v}_j(i - 1) + c_1r_1\left[\vec{P}_{\text{best},j} - \vec{x}_j(i-1)\right] + c_2r_2\left[\vec{G}_{\text{best},j} - \vec{x}_j(i - 1)\right]\]
  where
  \begin{align*}
    \vec{P}_{\text{best},j} &= \text{historical best position of } \vec{x}_j \\
    \vec{G}_{\text{best},j} &= \text{historical best position the entire swarm} \\
    c_1 &= \text{cognitive learning rate, typically 2} \\
    c_2 &= \text{social learning rate, typically 2} \\
    r_1,r_2 &= \text{random number} \in [0,1]
  \end{align*}
  then
  \[\vec{x}_j(i) = \vec{x}_j(i - 1) + \vec{v}_j(i)\]
\item If not converged, go to step~\ref{itm:swarm-update}.
\end{enumerate*}

\subsection*{Inertia}

This is an improvement to keep particle velocities from building up too fast.
\[\theta(i) = \theta_\mrm{max} - \left(\frac{\theta_\mrm{max} - \theta_\mrm{min}}{i_\mrm{max}}\right) i\]
where typically \(\theta_\mrm{min} = 0.4\) and \(\theta_\mrm{max} = 0.9\). Then, the velocity update relation is:
\[\vec{v}_j(i) = \theta(i - 1)\vec{v}_j(i - 1) + c_1r_1\left[\vec{P}_{\text{best},j} - \vec{x}_j(i-1)\right] + c_2r_2\left[\vec{G}_{\text{best},j} - \vec{x}_j(i - 1)\right]\]

\subsection*{Constraints}

This method works for inequality constraints only. Psuedo-\(F\):
\begin{description*}
\item[Stationary] fixed parameters, function of violation degree
\item[Nonstationary] dynamic penalty parameters
  \[F(\vec{x}) = f(\vec{x}) + C(i)H(\vec{x})\]
  where
  \begin{align*}
    C(i) &= {(ci)}^\alpha = \text{dynamic penalty parameter} \\
    H(\vec{x}) &= \sum_{j=1}^{m}{\left\{\psi\left[g_j(\vec{x})\right]{\left[q_j(\vec{x})\right]}^{\gamma\left[q_j(\vec{x})\right]}\right\}} = \text{penalty factor} \\
    \psi\left[g_j(\vec{x})\right] &= a\left(1 - \frac{1}{\mrm{e}^{q_j(\vec{x})}}\right) + b \\
    q_j(\vec{x}) &= \max\left(0, g_j(\vec{x})\right) \\
  \end{align*}
  where typically
  \begin{align*}
    c &= 0.5 \\
    \alpha &= 2 \\
    a &= 150 \\
    b &= 10 \\
    \gamma\left[g_j(\vec{x})\right] &= \begin{cases}
      1 & q_j(\vec{x}) \le 1 \\
      2 & q_j(\vec{x}) > 1
    \end{cases}
  \end{align*}
\end{description*}

\begin{thebibliography}{99}
\bibitem{flec} Ferguson, S. M., 2014, \emph{Engineering Design Optimization}, MAE 531 course
  lectures, North Carolina State University, Raleigh, NC.
\bibitem{eotap} Rao, S. S., 2009, \emph{Engineering Optimization: Theory and Practice}, 4th ed.,
  John Wiley \& Sons, Inc., Hoboken, NJ.
\end{thebibliography}

\end{document}