About Constraint Handling
An important difference between the available optimization solvers is the types of constraints they can handle, and how they do it. Constraints specified in the Optimization interface and Optimization study step can be divided into three categories:
•
Simple bounds are upper and lower bounds prescribed directly on the individual control variables, for example in the Optimization study step.
•
Pointwise constraints specify limits on an expression to be enforced at every node point in some region in space. In order to avoid excessively expensive gradient computations, such constraints are required to only depend on control variables directly and not indirectly via PDE solution variables. The constraint expression can, however, be a nonlinear expression in the control variables.
•
General constraints specify limits on global scalar expressions, typically evaluated as integrals over some domain. This generates a single constraint, as compared to one for each mesh node for a pointwise constraint. Therefore, the solvers can afford to compute a complete gradient also when the constraint is a, possible nonlinear, function of the PDE solution.
Note that all constraints are treated as inequalities. An equality constraint can be implemented by specifying the same upper and lower bound for an expression. However, not all constraint handling methods are able to deal with the reduction in control variable space dimension which this implies. Therefore, when possible, it is better to perform a manual change of variables, eliminating a control variable dimension.
Constraint Handling for Derivative-Free Methods
The derivative-free methods Nelder–Mead, COBYLA, and EGO can internally handle simple bounds and general constraints on global scalar expressions. In practice this means that in addition to simple bounds and constraints defined in the Optimization study step, Integral Inequality Constraint nodes and Global Inequality Constraint nodes in Optimization interfaces are accounted for.
The constraint handling algorithm used in Nelder–Mead is in principle based on filtering out candidate points in the control variable space which fall outside the feasible region, and to some extent adjust search directions accordingly. The procedure is not guaranteed to find a constrained local minimum fulfilling the KKT conditions.
COBYLA, in contrast, approximates objective function and constraints in a uniform way. Therefore, provided all functions are sufficiently smooth, it will in general find an approximate constrained local minimum. The returned solution may, however, lie slightly outside the feasible set. This can happen, in particular, if the constraints are nonlinear.
The constraint handling algorithm in EGO is based on incorporating the constraints directly into the surrogate model using Gaussian process functions. The algorithm models both the objective function and constraints, balancing their improvement when selecting new sample points using the acquisition function. While the procedure improves constraint satisfaction as the model evolves, the returned solution may still lie slightly outside the feasible region, especially in the early stages or if the constraints are highly nonlinear.
The Augmented Lagrangian Method
BOBYQA handles simple bounds internally, but general constraints only via an external iterative procedure based on repeated minimization of an augmented Lagrangian. This augmented Lagrangian method can also be used as an alternative to the internal penalty method in the Nelder–Mead solver, but is not selected by default.
The basic principle behind the augmented Lagrangian method is to include the Lagrange multipliers of general constraints as control variables in an augmented problem. In the first iteration, the Lagrange multipliers are set to zero and a modified objective function including a quadratic penalty for constraint violation is minimized. This gives a solution which is in general outside the feasible set, that is, it violates the constraints. In each subsequent iteration, the Lagrange multipliers — as well as a number of penalty parameters — are updated based on the current constraint violation, and a new subproblem is solved. The sequence of subproblems, which converges toward the feasible set from the outside, is terminated once a specified tolerance for constrain violation is reached.
The augmented Lagrangian method as such is very general and quite robust, but to be efficient, it requires balancing the effort spent on each subproblem against the improvement in each outer iteration. It also requires selection of an initial penalty factor. A higher penalty factor generally leads to faster convergence of the augmented Lagrangian algorithm, but subproblems become more ill-conditioned and there is a threshold beyond which the method may become unstable; conversely, a lower penalty factor makes the algorithm more robust but required more iterations.
In the Settings window for the Optimization study step, you can choose to use an Automatic or a Manual definition of the initial Penalty parameter ρ. The automatic setting is default and computes an initial value based on the objective and constraint function values at the initial point. You can also limit the Maximum number of augmented iterations and select a strategy for updating δ (tolerance for the subsolver). The options Dynamic I and Dynamic II both tighten the subsolver tolerance from iteration to iteration, the latter providing some additional control. There is also a Manual option. Finally, specify the Constraint tolerance, that is, the maximum allowable constraint violation in the final solution.
Since the augmented Lagrangian method computes Lagrange multipliers for each constraint explicitly, these are also available for postprocessing. Their values represent the sensitivity (derivative) of the objective function with respect to changes in a constraint bound. In the Insert Expression and Add Expression menus, available for most postprocessing features, you will find the Lagrange multipliers under Model > Solver > Lagrange multipliers.
Constraint Handling for Gradient-Based Methods
The SNOPT algorithm handles constraints of all types efficiently. Constraint handling in this SQP method is based on linearizing the constraint in the outer, major, iteration, and using an active-set QP solver in the inner loop to decide which constraints are active and bounding the solution at the current iterate. This process requires accurate evaluation of the gradient of the constraints, also known as the constraint Jacobian.
The IPOPT algorithm supports constraints and bounds by using an interior point line search method, see IPOPT-Specific Settings.
The MMA algorithm accepts constraints of the same general type as SNOPT, requiring an accurate constraint Jacobian, but handles them differently. In each outer, major, iteration, linear and linearized constraints are combined with a linearized objective function into a convex smooth approximation whose unique optimum is always feasible unless the feasible set is empty. The globally convergent version of MMA implemented in the Optimization module is conservative in a way which ensures that each major iterate is feasible not only with respect to the linearized constraints, but with respect to the fully nonlinear constraints.
In the MMA implementation, simple bounds are much less expensive memory-wise than global or pointwise constraints. In particular for control variable fields, simple control variable bounds (added by right-clicking on the Control Variable Field feature) are more efficient than enforcing the same bounds using a Pointwise Inequality Constraint feature.
The Levenberg–Marquardt solver supports bounds but not constraints.