2.1. Review of Optimality Criteria Method
In this section, the conventional OCM [
1] is reviewed for the purpose of developing the GOCM in the following section. In the SIMP method, the topological density of each element,
, is considered as a design variable. Then, the optimization problem can be stated as
In Equation (1),
is the vector of design variables,
is the compliance,
is the penalization power (typically
) in the SIMP method,
is the material volume with
being the volume of the element,
is the design domain volume, and
is the volume fraction.
is the stiffness matrix of element
e and nodal degrees of freedom (DOFs)
. The assembly of the element stiffness matrix yields the global stiffness matrix
, and that of nodal DOFs yields the global vector of DOFs
[
18]. The topological densities have the upper and lower bounds,
and
, respectively. Many topology optimization problems are defined in the rectangular grid of mesh, i.e., pixels in 2D and voxels in 3D. However, the optimization problem in Equation (1) can also be applied to irregular grids.
In Equation (1), since the structural equilibrium
is solved first for given design variables, it is automatically satisfied at each design iteration. Moreover, the side constraints,
, can be satisfied directly when design variables are determined. Therefore, the optimization problem has a single compliance objective function and a single constraint of the volume fraction. The particular optimization problem is attractive because the sensitivity comes almost free of computation. The sensitivities of the compliance and the volume fraction can respectively be calculated as
The sensitivity of the compliance is so-called self-adjoint, which means that the adjoint response is identical to the structural response. Therefore, no additional calculation is required for the adjoint response. The sensitivity of the volume fraction is nothing but the element volume itself.
The OCM can be derived by converting the constrained optimization problem in Equation (1) into an unconstrained one by defining the following Lagrange function:
The Karush–Kuhn–Tucker first-order optimality condition becomes
The procedure of OCM in Bendsøe [
1] is composed of two-level loops. In the inner loop, the design variable
is updated to satisfy the first condition in Equation (5) for a given Lagrange multiplier
. In the outer loop, the Lagrange multiplier is updated to satisfy the volume fraction constraint. More specifically, the OCM changes the design variable in such a way that the first equation in Equation (5) becomes zero. For that purpose, the following scale factor is defined for each element:
Since the optimality condition is satisfied when the objective sensitivity is equal and opposite of constraint sensitivity multiplied by the Lagrange multiplier [
15], the first equation in the optimality condition is satisfied when
. The idea of OCM is to change design variables based on the scale factor as
The design would not be changed when
because the optimality condition is already satisfied. When
, it means that increasing design
is less efficient in decreasing the compliance than increasing the volume. In this case, therefore, it is better to reduce
. When
, the opposite is true, and the design variable should increase. In addition, since it is not preferred to change a design significantly in one iteration, the maximum change in design
is also imposed.
Figure 1 illustrates the possible ranges of design change. Starting from the current design
, the first case is when
. In this case, the design can be changed within the maximum change. The second case is when
. In this case, the design can be changed in
. When
, the design can be changed in
.
In the outer loop, the Lagrange multiplier is determined to satisfy the volume fraction constraint. The Matlab implementation of the 99-line topology optimization code [
13] used a bisection method to find the Lagrange multiplier. Starting from the lower and upper bounds [0, 100,000] of the Lagrange multiplier, the range is halved every iteration of the outer loop, and the Lagrange multiplier takes the value in the middle of the range. With the current Lagrange multiplier, if the constraint is positive,
, then the upper half of the range is used in the next iteration; otherwise, the lower half is used. This bisection process is repeated until the range becomes less than a convergence tolerance.
For the purpose of generalizing the OCM, three obstacles must be overcome. The first obstacle is the compliance objective and the volume fraction constraint. First, the sensitivity of volume in Equation (3) is constant and independent of design. In addition, the sensitivity of compliance in Equation (2) can be calculated once the structural equilibrium is solved for . Therefore, this particular combination of the objective and constraint has the almost free computation of sensitivities. However, the procedure cannot be generalized when there is more than one constraint.
The second obstacle is the assumption of sensitivities. In Equation (6), the OCM algorithm assumes that the objective derivative is negative and the constraint derivative is positive. This is true with the compliance objective and the volume fraction constraint. From Equation (2), the quadratic form of a positive semi-definite stiffness matrix is always non-negative. Therefore, the compliance sensitivity in Equation (2) is always negative or zero. On the other hand, the volume fraction sensitivity in Equation (3) is constant and positive. In general optimization problems, it is possible that the objective may have a positive sensitivity for some designs, while a negative sensitivity for others. In particular, when there exists a design-dependent load, such as gravity, the sensitivity of compliance can be negative for some elements. Therefore, the design update formula in Equation (7) cannot be used as it is.
The third obstacle is the computational cost related to finding the Lagrange multiplier. As mentioned before, the Lagrange multiplier is determined through the bisection method, which requires hundreds of outer-loop iterations. When the number of design variables is large, this process can take a lot of computational time. The same is true for the method of moving asymptotes (MMA), where solving the MMA subproblem can be expensive especially when multiple constraints are active.
2.2. Generalized Optimality Criteria Method
The generalized optimality criteria (GOCM) method for topology optimization extends the capability of the OCM to multiple inequality constraints, possibly with improved computational efficiency. The general idea of OCM has been available for a long time in parameter optimization, albeit it was limited to minimizing weight with displacements, stresses, and natural frequencies [
16,
17]. The basic idea is to solve the Karush–Kuhn–Tucker condition, which is the necessary condition for optimization. The general optimization problem can be defined as
Even if a single objective function and only less-than-or-equal-to-type inequality constraints are used, it can easily be generalized to multiple objective functions with weighted sum and other types of constraints. In the following derivations, it is assumed that both the objective and constraints are normalized using the initial value and constraint bounds. For example, in the case of stress constraint given as can be normalized as . In the case of the objective function, it can be normalized using the initial value.
In general, the constrained optimization problem can be converted into an unconstrained optimization problem using either the Lagrange multiplier method or penalty method. In the Lagrange multiplier method, the Lagrange function is defined by combining the objective function and constraints using Lagrange multipliers as
where
is the Lagrange multiplier corresponding to constraint
.
is called a slack variable, which is not zero when the constraint is inactive (i.e., less than zero).
The necessary condition for optimum is when the Lagrange function is stationary, i.e., its derivatives are zero. Since the Lagrange function has three variables, it is differentiated by all three variables as
where
is the column vector of gradients. Due to the complementary slackness (i.e., switching condition),
, only the active constraints need to be considered in the necessary condition.
Since the second and third equations in Equation (10) are satisfied by identifying active constraints, the process of GOCM is to solve the first part of Equation (10). At an optimum design, the objective sensitivity can be represented by a linear combination of active constraint gradients. The coefficients in the linear combination are indeed the Lagrange multipliers. Since the Lagrange multipliers,
, and design variables,
, are coupled, they have to be solved simultaneously. The challenge is that solving the nonlinear equation is computationally intensive and difficult due to numerical instability. The initial implementation of OCM by Sigmund [
13] has a double-loop method, where the inner loop calculates the design variable, while the outer loop calculates the Lagrange multiplier. When multiple constraints are active, however, it would require multiple levels of loops to calculate the coupled equations.
The novel idea of the proposed GOCM is that the Lagrange multipliers do not have to satisfy Equation (10) at every iteration. Therefore, no iterative bisection is required to find the Lagrange multipliers. Instead, the Lagrange multipliers are updated at every optimization step so that when the optimization is converged, the necessary condition in Equation (10) is satisfied. Patnaik et al. [
16] proposed several methods of updating the Lagrange multiplier.
Table 1 summarizes the three updating schemes. In the table,
is the initial value of an update parameter, and
is an acceleration parameter used to modify the update parameter. Patnaik et al. [
16] suggested
and
.
In this article, the updating algorithm of the Lagrange multiplier is composed of two steps: (a) initial estimate and (b) update during iteration. First, the initial values of Lagrange multipliers are estimated using the inverse form in
Table 1. However, it would be computationally complicated to calculate the inverse of the constraint gradient matrix, and it would be unnecessary as it provides the initial estimate. Instead, the uncoupled version of the inverse form is used for the initial estimate as
The above estimate works well in most cases except when either the objective gradient or constraint gradients are zero or too small. In the optimization problem formulation, the objective function and constraints are normalized to the initial value and/or constraint limit. Therefore, the initial value of is a good starting point when .
The basic concept of GOCM is when the constraint is violated, the Lagrange multiplier needs to increase, while it is decreased when the constraint is inactive. Therefore, a simple updating formula would be
. In this simple formula, the Lagrange multiplier increases when the constraint is positive (i.e., the constraint is violated) or it decreases when the constraint is negative (i.e., the constraint is inactive). This simple algorithm can cause some problems when the constraint violation is not recovered quickly, or the constraint is inactive for many steps. For example, if the constraint is violated a lot and stays violated for many steps, then the Lagrange multipliers can increase too fast. Therefore, it would be necessary to include the effect of change in constraint. In the article, the following updating rule is proposed for the Lagrange multiplier:
where
is the update parameter similar to those used in
Table 1. The update parameter can have a positive value when both the sign of
and
are the same. A more complex scheme can also be used, as shown in the Matlab code in
Section 3.
The updating algorithm in Equation (12) is applied (a) when the current constraint is violated and the constraint increases from the previous step or (b) when the current constraint is inactive and the constraint continues to decrease. In order to make the updating process stable, the maximum change of the Lagrange multipliers is used. Moreover, the Lagrange multipliers are limited to change within the lower- and upper-bounds. Since all constraints are normalized by their bounds, the magnitudes of the Lagrange multipliers are in a similar order of magnitude.
In the theory of Lagrange multiplier, it is well known that it is zero when the constraint is not active, while it is positive when the constraint is active. That means, it is unnecessary to consider the Lagrange multiplier for inactive constraints and only active constraints are considered in the optimization. However, from a practical point of view, it is difficult to turn on and turn off constraints during optimization. Therefore, in the implementation, all constraints and Lagrange multipliers are retained. Then a constraint is inactive, the corresponding Lagrange multiplier will converge to its lower bound, which reduces its effect on the optimality criteria.
Once the Lagrange multipliers are updated, the next step of GOCM is to update design variables. When the optimization problem is composed of the compliance objective and volume fraction constraint, the scale factor
in Equation (6) is used to update design variables. This updating algorithm is based on the assumption that that the objective sensitivity is negative, while the constraint sensitivity is positive. Similar algorithms for design variable update were available in Patnaik et al. [
16], which are summarized in
Table 2. All the formulas are based on the scale factor defined as
Making
is equivalent to the stationary condition of Karush–Kuhn–Tucker condition in Equation (10). The two acceleration parameters in
Table 2 are suggested to be
and
.
In order to consider the general objective and constraints, we need to modify the scale factor in Equation (13). For a given design variable (i.e., element), the scale factor is calculated on the basis of the sign of sensitivities. The numerator has all terms with negative sensitivities, while the denominator has all terms with positive sensitivities. Accordingly, the scale factor is calculated as
where
and
. When
, this formula will also satisfy the stationary condition of the Lagrange function. The only difference is that the original formulation calculates the scale factor on the basis of the ratio between objective sensitivity and the weighted sum of constraint sensitivities, while the proposed method uses the ratio between the positive and negative sensitivities. In the case of compliance objective and volume fraction constraint, Equations (6), (13), and (14) yield the identical scale factor.
In some special situations, the scale factor needs to be modified. If the numerator and/or the denominator are zero, the scale factor needs to be limited so that it stays close to one. Moreover, is scaled such that the design change in each iteration is less than . Once the scale factor is determined, Equation (7) is used to update the design variable.