Author

Philip E. Gill; Department of Mathematics, University of California, San Diego, La Jolla, CA

Walter Murray, Michael A. Saunders; Department of EESOR, Stanford University, Stanford, CA

Introduction

This section describes the GAMS interface to the general-purpose NLP solver SNOPT, (Sparse Nonlinear Optimizer) which implements a sequential quadratic programming (SQP) method for solving constrained optimization problems with smooth nonlinear functions in the objective and constraints. The optimization problem is assumed to be stated in the form

\begin{equation} \tag{1} \begin{array}{rllll} \textrm{NP}: & ~~~~~~~~~~~~ & \textrm{minimize} \> \textrm{or} \> \textrm{maximize} \> f_0(x) & & \\[5pt] & ~~~~~~~~~~~~ & \textrm{subject to} \> {\begin{array}{l} f(x) \sim b_1\\A_L x \sim b_2 \\l \le x \le u, \end{array}} & & \end{array} \end{equation}

where \(x \in \Re^n\), \(f_0(x)\) is a linear or nonlinear smooth objective function, \(l\) and \(u\) are constant lower and upper bounds, \(f(x)\) is a set of nonlinear constraint functions, \(A_L\) is a sparse matrix, \(\sim\) is a vector of relational operators ( \(\le\), \(\ge\) or \(=\)), and \(b_1\) and \(b_2\) are right-hand side constants. \(f(x)\sim b_1\)s are the nonlinear constraints of the model and \(A_L x \sim b_2\) form the linear constraints.

The gradients of \(f_0\) and \(f_i\) are automatically provided by GAMS, using its automatic differentiation engine.

The bounds may include special values -INF or +INF to indicate \(l_j = -\infty\) or \(u_j=+\infty\) for appropriate \(j\). Free variables have both bounds infinite and fixed variables have \(l_j = u_j\).

SNOPT uses a sequential quadratic programming (SQP) algorithm that obtains search directions from a sequence of quadratic programming subproblems. Each QP subproblem minimizes a quadratic model of a certain Lagrangian function subject to a linearization of the constraints. An augmented Lagrangian merit function is reduced along each search direction to ensure convergence from any starting point.

SNOPT is most efficient if only some of the variables enter nonlinearly, or if the number of active constraints (including simple bounds) is nearly as large as the number of variables. SNOPT requires relatively few evaluations of the problem functions.

Problem Types

If the nonlinear functions are absent, the problem is a linear program (LP) and SNOPT applies the primal simplex method [40] . Sparse basis factors are maintained by LUSOL [81] as in MINOS [140] .

If only the objective is nonlinear, the problem is linearly constrained (LC) and tends to solve more easily than the general case with nonlinear constraints (NC). Note that GAMS models have an objective variable instead of an objective function. The GAMS/SNOPT link will try to substitute out the objective variable and reformulate the model such that SNOPT will see a true objective function.

For both linearly and nonlinearly constrained problems SNOPT applies a sparse sequential quadratic programming (SQP) method [84] using limited-memory quasi-Newton approximations to the Hessian of the Lagrangian. The merit function for steplength control is an augmented Lagrangian, as in the dense SQP solver NPSOL [80, 83] .

In general, SNOPT requires less matrix computation than NPSOL and fewer evaluations of the functions than the nonlinear algorithms in MINOS [138, 139] . It is suitable for nonlinear problems with thousands of constraints and variables, and is most efficient if only some of the variables enter nonlinearly, or there are relatively few degrees of freedom at a solution (i.e., many constraints are active). However, unlike previous versions of SNOPT, there is no limit on the number of degrees of freedom.

Selecting the SNOPT Solver

If SNOPT is not specified as the default solver for the desired model type (e.g. NLP), then the following statement can be used in your GAMS model to select SNOPT:

option nlp=SNOPT;

The option statement should appear before the solve statement. To be complete, we mention that the solver can be also specified on the command line, as in:

> gams camcge nlp=snopt

This will override the global default, but if an algorithm option has been specified inside the model, then that specification takes precedence.

If the model contains non-smooth functions like abs \((x)\), or \(\max(x,y)\) you can try to get it solved by SNOPT using

option dnlp=SNOPT;

These models have discontinuous derivatives however, and SNOPT was not designed for solving such models. Discontinuities in the gradients can sometimes be tolerated if they appear away from an optimum.

Description of the method

Here we briefly describe the main features of the SQP algorithm used in SNOPT and introduce some terminology. The SQP algorithm is fully described by by Gill, Murray and Saunders [85] .

Objective function reconstruction

The first step GAMS/SNOPT performs is to try to reconstruct the objective function. In GAMS, optimization models minimize or maximize an objective variable. SNOPT however works with an objective function. One way of dealing with this is to add a dummy linear function with just the objective variable. Consider the following GAMS fragment:

    obj.. z =e= sum(i, sqr[r(i)]);

    model m /all/;
    solve m using nlp minimizing z;

This can be cast in form (1) by saying minimize \(z\) subject to \(z = \sum_i r^2_i\) and the other constraints in the model. Although simple, this approach is not always preferable. Especially when all constraints are linear it is important to minimize the nonlinear expression \(\sum_i r^2_i\) directly. This can be achieved by a simple reformulation: \(z\) can be substituted out. The substitution mechanism carries out the formulation if all of the following conditions hold:

• the objective variable \(z\) is a free continuous variable (no bounds are defined on \(z\)),
• \(z\) appears linearly in the objective function,
• the objective function is formulated as an equality constraint,
• \(z\) is only present in the objective function and not in other constraints.

For many models it is very important that the nonlinear objective function be used by SNOPT. For instance the model chem.gms from the model library solves in 16 iterations. When we add the bound

   energy.lo = 0;

on the objective variable energy and thus preventing it from being substituted out, SNOPT will not be able to find a feasible point for the given starting point.

This reformulation mechanism has been extended for substitutions along the diagonal. For example, the GAMS model

   variables x,y,z;
   equations e1,e2;
   e1..z =e= y;
   e2..y =e= sqr(1+x);
   model m /all/;
   option nlp=snopt;
   solve m using nlp minimizing z;

will be reformulated as an unconstrained optimization problem

\[ \min f(x) = (1+x)^2. \]

Constraints and slack variables

Problem (1) contains \(n\) variables in \(x\). Let \(m\) be the number of components of \(f(x)\) and \(A_L x\) combined. The upper and lower bounds on those terms define the general constraints of the problem. SNOPT converts the general constraints to equalities by introducing a set of slack variables \(s = (s_1, s_2, ..., s_m)^T\). For example, the linear constraint \(5 \le 2x_1 + 3x_2 \le +\infty\) is replaced by \(2x_1 + 3x_2 - s_1 = 0\) together with the bounded slack \(5 \le s_1 \le +\infty\). Problem (1) can be written in the equivalent form

\[ \begin{array}{l@{\hspace{1ex}}l} \textrm{minimize} & \ \ \ f_0(x) \\ \textrm{subject to} & {\left( \begin{matrix} f(x) \\ A_L x \end{matrix}\ \right)} - s =0, \quad l \le {\left( \begin{matrix} x \\ s \end{matrix}\ \right)} \le u. \end{array} \]

where a maximization problem is cast into a minimization by multiplying the objective function by \(-1\).

The linear and nonlinear general constraints become equalities of the form \(f(x) - s_N = 0\) and \(A_L x - s_L = 0\), where \(s_L\) and \(s_N\) are known as the linear and nonlinear slacks.

Major iterations

The basic structure of SNOPT's solution algorithm involves major and minor iterations. The major iterations generate a sequence of iterates \((x_k)\) that satisfy the linear constraints and converge to a point that satisfies the first-order conditions for optimality. At each iterate \(\{x_k\}\) a QP subproblem is used to generate a search direction towards the next iterate \(\{x_{k+1}\}\). The constraints of the subproblem are formed from the linear constraints \(A_L x - s_L =0\) and the nonlinear constraint linearization

\[ f(x_k) + f'(x_k)(x - x_k) - s_N = 0, \]

where \(f'(x_k)\) denotes the Jacobian: a matrix whose rows are the first derivatives of \(f(x)\) evaluated at \(x_k\). The QP constraints therefore comprise the \(m\) linear constraints

\[ \begin{array}{rrr@{\hspace{4pt}}c@{\hspace{4pt}}l} f'(x_k)x & - s_N & &=& -f(x_k) + f'(x_k) x_k, \\[1ex] A_L x & & -s_L &=& \> 0, \end{array} \]

where \(x\) and \(s\) are bounded by \(l\) and \(u\) as before. If the \(m\times n\) matrix \(A\) and \(m\)-vector \(b\) are defined as

\[ A = {\left(\! \begin{matrix} f'(x_k) \\ A_L \end{matrix}\! \right)} ~~\text{and}~~ b = {\left(\! \begin{matrix} -f(x_k) + f'(x_k) x_k \\ 0 \end{matrix}\! \right)}, \]

then the QP subproblem can be written as

\begin{equation} \tag{2} \begin{array}{rllll} \textrm{QP}_k : & ~~~~~~~~~~~~ & \displaystyle\min_{{x,s}} & {\displaystyle q(x,x_k)=g_k^T(x-x_k)+\frac{1}{2}(x-x_k)^T H_k (x-x_k)} & \\[5pt] & ~~~~~~~~~~~~ & \textrm{subject to} & {\displaystyle Ax-s=b, \ \ l \le {\left( \begin{matrix} x \\ s \end{matrix}\ \right)} \le u, } & \end{array} \end{equation}

where \(q(x,x_k)\) is a quadratic approximation to a modified Lagrangian function [84] . The matrix \(H_k\) is a quasi-Newton approximation to the Hessian of the Lagrangian. A BFGS update is applied after each major iteration. If some of the variables enter the Lagrangian linearly the Hessian will have some zero rows and columns. If the nonlinear variables appear first, then only the leading \(n_1\) rows and columns of the Hessian need be approximated, where \(n_1\) is the number of nonlinear variables.

Minor iterations

Solving the QP subproblem is itself an iterative procedure. Here, the iterations of the QP solver SQOPT [86] form the minor iterations of the SQP method.

SQOPT uses a reduced-Hessian active-set method implemented as a reduced-gradient method similar to that in MINOS [138] .

At each minor iteration, the constraints \(A x - s = b\) are partitioned into the form

\[ Bx_B + Sx_S + N x_N = b, \]

where the basis matrix \(B\) is square and nonsingular and the matrices \(S\), \(N\) are the remaining columns of \((A\> {-}I)\). The vectors \(x_B\), \(x_S\), \(x_N\) are the associated basic, superbasic, and nonbasic components of the variables \((x, s)\).

The term active-set method arises because the nonbasic variables \(x_N\) are temporarily frozen at their upper or lower bounds, and their bounds are considered to be active. Since the general constraints are satisfied also, the set of active constraints takes the form

\[ {\left(\begin{matrix} B & S & N \\ & & I \end{matrix} \right)} {\left(\begin{matrix} x_B \\ x_S \\ x_N \end{matrix} \right)} = {\left(\begin{matrix} b \\ x_N \end{matrix} \right)}, \]

where \(x_N\) represents the current values of the nonbasic variables. (In practice, nonbasic variables are sometimes frozen at values strictly between their bounds.) The reduced-gradient method chooses to move the superbasic variables in a direction that will improve the objective function. The basic variables "tag along" to keep \(Ax-s = b\) satisfied, and the nonbasic variables remain unaltered until one of them is chosen to become superbasic.

At a nonoptimal feasible point \((x, s)\) we seek a search direction \(p\) such that \((x, s) + p\) remains on the set of active constraints yet improves the QP objective. If the new point is to be feasible, we must have \(Bp_B + Sp_S + Np_N = 0\) and \(p_N = 0\). Once \(p_S\) is specified, \(p_B\) is uniquely determined from the system \(Bp_B = -Sp_S\). It follows that the superbasic variables may be regarded as independent variables that are free to move in any desired direction. The number of superbasic variables ( \(n_S\) say) therefore indicates the number of degrees of freedom remaining after the constraints have been satisfied. In broad terms, \(n_S\) is a measure of how nonlinear the problem is. In particular, \(n_S\) need not be more than one for linear problems.

The reduced Hessian and reduced gradient

The dependence of \(p\) on \(p_S\) may be expressed compactly as \(p = Zp_S\), where \(Z\) is a matrix that spans the null space of the active constraints:

\begin{equation} \tag{3} Z = P {\left( \begin{matrix} -B^{-1}S \\ I \\ 0 \end{matrix} \right)} \end{equation}

where \(P\) permutes the columns of \((A\> {-}I)\) into the order \((B\> S\> N)\). Minimizing \(q(x, x_k)\) with respect to \(p_S\) now involves a quadratic function of \(p_S\):

\begin{equation} \tag{4} g^TZp_S + \frac{1}{2}p^T_S Z^T H Z p_S, \end{equation}

where \(g\) and \(H\) are expanded forms of \(g_k\) and \(H_k\) defined for all variables \((x, s)\). This is a quadratic with Hessian \(Z^THZ\) (the reduced Hessian) and constant vector \(Z^Tg\) (the reduced gradient). If the reduced Hessian is nonsingular, \(p_S\) is computed from the system

\begin{equation} \tag{5} \label{redhessian} Z^THZp_S = -Z^Tg. \end{equation}

The matrix \(Z\) is used only as an operator, i.e., it is not stored explicitly. Products of the form \(Zv\) or \(Z^Tg\) are obtained by solving with \(B\) or \(B^T\) . The package LUSOL [81] is used to maintain sparse \(LU\) factors of \(B\) as the \(BSN\) partition changes. From the definition of \(Z\), we see that the reduced gradient can be computed from

\[ B^T\pi = g_B,\>\> Z^Tg = g_S - S^T\pi, \]

where \(\pi\) is an estimate of the dual variables associated with the \(m\) equality constraints \(Ax - s = b\), and \(g_B\) is the basic part of \(g\).

By analogy with the elements of \(Z^Tg\), we define a vector of reduced gradients (or reduced costs) for all variables in \((x, s)\):

\[ d = g - {\left(\begin{matrix} A^T\\ {-}I \end{matrix} \right)} \pi, \ \ \text{so that $d_S = Z^Tg$.} \]

At a feasible point, the reduced gradients for the slacks \(s\) are the dual variables \(\pi\).

The optimality conditions for subproblem QP \(_k\) (2) may be written in terms of \(d\). The current point is optimal if \(d_j \ge 0\) for all nonbasic variables at their lower bounds, \(d_j \le 0\) for all nonbasic variables at their upper bounds, and \(d_j = 0\) for all superbasic variables ( \(d_S = 0\)). In practice, SNOPT requests an approximate QP solution \(({\skew{2.8}\widehat x}_k, {\widehat s}_k, {\skew1\widehat \pi}_k)\) with slightly relaxed conditions on \(d_j\).

If \(d_S = 0\), no improvement can be made with the current \(BSN\) partition, and a nonbasic variable with non-optimal reduced gradient is selected to be added to \(S\). The iteration is then repeated with \(n_S\) increased by one. At all stages, if the step \((x, s) + p\) would cause a basic or superbasic variable to violate one of its bounds, a shorter step \((x, s) + \alpha p\) is taken, one of the variables is made nonbasic, and \(n_S\) is decreased by one.

The process of computing and testing reduced gradients \(d_N\) is known as pricing (a term introduced in the context of the simplex method for linear programming). Pricing the \(j\)th variable means computing \(d_j = g_j - a^T_j \pi\), where \(a_j\) is the \(j\)th column of \((A\> {-}I)\). If there are significantly more variables than general constraints (i.e., \(n \gg m\)), pricing can be computationally expensive. In this case, a strategy known as partial pricing can be used to compute and test only a subset of \(d_N\).

Solving the reduced Hessian system (5) is sometimes expensive. With the option QPSolver Cholesky, an upper-triangular matrix \(R\) is maintained satisfying \(R^TR = Z^THZ\). Normally, \(R\) is computed from \(Z^THZ\) at the start of phase 2 and is then updated as the \(BSN\) sets change. For efficiency the dimension of \(R\) should not be excessive (say, \(n_S \le 1000\)). This is guaranteed if the number of nonlinear variables is "moderate". Other QPSolver options are available for problems with many degrees of freedom.

The merit function

After a QP subproblem has been solved, new estimates of the NLP solution are computed using a linesearch on the augmented Lagrangian merit function

\begin{equation} \tag{6} \label{eqn-def-merit} \mathcal{M}(x,s,\pi) = f(x) - \pi^T \bigl( F(x) - s_N \bigr) + \frac{1}{2} \bigl( F(x) - s_N \bigr)^T D \bigl( F(x) - s_N\bigr), \end{equation}

where \(D\) is a diagonal matrix of penalty parameters. If \((x_k,s_k,\pi_k)\) denotes the current solution estimate and \((\skew{2.8}\widehat x_k, \widehat s_k, \skew1\widehat \pi_k)\) denotes the optimal QP solution, the linesearch determines a step \(\alpha_k\) ( \(0 < \alpha_k \le 1\)) such that the new point

\begin{equation} \tag{7} \label{eqn-def-newvars} {\left(\begin{matrix} x_{k+1} \\ s_{k+1} \\ \pi_{k+1} \end{matrix} \right)} = {\left(\begin{matrix} x_k \\ s_k \\ \pi_k \end{matrix} \right)} + \alpha_k {\left(\begin{matrix} \skew{2.8} \widehat x_k - x_k \\ \widehat s_k - s_k \\ \skew1 \widehat \pi_k - \pi_k \end{matrix} \right)} \end{equation}

gives a sufficient decrease in the merit function. When necessary, the penalties in \(D\) are increased by the minimum-norm perturbation that ensures descent for \(\mathcal{M}\) [83] . As in NPSOL, \(s_N\) is adjusted to minimize the merit function as a function of \(s\) prior to the solution of the QP subproblem. For more details, see [80, 53] .

Treatment of constraint infeasibilities

SNOPT makes explicit allowance for infeasible constraints. Infeasible linear constraints are detected first by solving a problem of the form

\[ \begin{array}{rllll} \textrm{FLP}: & ~~~~~~~~~~~~ & \textrm{minimize} & e^T(v + w) & \\[5pt] & ~~~~~~~~~~~~ & \textrm{subject to} & \ell \le \left( \begin{matrix} x\\ A_L x-v+w \end{matrix} \right) \le u, \> v \ge 0,\> w \ge 0, & \end{array} \]

where \(e\) is a vector of ones. This is equivalent to minimizing the sum of the general linear constraint violations subject to the simple bounds. (In the linear programming literature, the approach is often called elastic programming. We also describe it as minimizing the \(\ell_1\) norm of the infeasibilities.)

If the linear constraints are infeasible ( \(v \ne 0\) or \(w \ne 0\)), SNOPT terminates without computing the nonlinear functions.

If the linear constraints are feasible, all subsequent iterates satisfy the linear constraints. (Such a strategy allows linear constraints to be used to define a region in which the functions can be safely evaluated.) SNOPT proceeds to solve NP (1) as given, using search directions obtained from a sequence of quadratic programming subproblems (2).

If a QP subproblem proves to be infeasible or unbounded (or if the dual variables \(\pi\) for the nonlinear constraints become large), SNOPT enters "elastic" mode and solves the problem

\[ \begin{array}{rllll} \textrm{NP}(\gamma): & ~~~~~~~~~~~~ & \textrm{minimize} & {f_0(x) + \gamma e^T(v + w)} & \\[5pt] & ~~~~~~~~~~~~ & \textrm{subject to} & {\ell \le \left( \begin{matrix} x\\ f(x)-v+w\\ A_Lx \end{matrix} \right) \le u, \> v \ge 0,\> w \ge 0,} & \end{array} \]

where \(\gamma\) is a nonnegative parameter (the elastic weight), and \(f(x) + \gamma e^T(v + w)\) is called a composite objective. If \(\gamma\) is sufficiently large, this is equivalent to minimizing the sum of the nonlinear constraint violations subject to the linear constraints and bounds. A similar \(\ell_1\) formulation of NP is fundamental to the S \(\ell_1\)QP algorithm of Fletcher [63] . See also Conn [37] .

The initial value of \(\gamma\) is controlled by the optional parameter elastic weight.

Starting points and advanced bases

A good starting point may be essential for solving nonlinear models. We show how such a starting point can be specified in a GAMS environment, and how SNOPT will use this information.

A related issue is the use of "restart" information in case a number of related models are solved in a row. Starting from an optimal point of a previous solve statement is in such situations often beneficial. In a GAMS environment this means reusing primal and dual information, which is stored in the .L and .M fields of variables and equations.

Starting points

To specify a starting point for SNOPT use the .L level values in GAMS. For example, to set all variables \(x_{i,j} := 1\) use x.l(i,j)=1;. The default values for level values are zero.

Setting a good starting point can be crucial for getting good results. As an (artificial) example consider the problem where we want to find the smallest circle that contains a number of points \((x_i,y_i)\):

\[ \begin{array}{rllll} \textrm{Example}: & ~~~~~~~~~~~~ & \textrm{minimize} & r & \\[5pt] & ~~~~~~~~~~~~ & \textrm{subject to} & (x_i-a)^2 + (y_i-b)^2 \le r^2, \> r \ge 0. & \end{array} \]

This problem can be modeled in GAMS as follows.

set i points /p1*p10/;

parameters
    x(i)    x coordinates,
    y(i)    y coordinates;

* fill with random data
x(i) = uniform(1,10);
y(i) = uniform(1,10);

variables
    a       x coordinate of center of circle
    b       y coordinate of center of circle
    r       radius;

equations
    e(i)    points must be inside circle;


e(i).. sqr(x(i)-a) + sqr(y(i)-b) =l= sqr(r);

r.lo = 0;

model m /all/;
option nlp=snopt;
solve m using nlp minimizing r;

Without help, SNOPT will not be able to find an optimal solution. The problem will be declared infeasible. In this case, providing a good starting point is very easy. If we define

\begin{eqnarray*} x_{\min} &=& \min_i x_i,\\ y_{\min} &=& \min_i y_i,\\ x_{\max} &=& \max_i x_i,\\ y_{\max} &=& \max_i y_i, \end{eqnarray*}

then good estimates are

\begin{eqnarray*} a &=& (x_{\min}+x_{\max})/2, \\ b &=& (y_{\min}+y_{\max})/2, \\ r &=& \sqrt{ (a-x_{\min})^2 + (b-y_{\min})^2}. \end{eqnarray*}

Thus we include in our model:

parameters xmin,ymin,xmax,ymax;
xmin = smin(i, x(i));
ymin = smin(i, x(i));
xmax = smax(i, x(i));
ymax = smax(i, y(i));

* set starting point
a.l = (xmin+xmax)/2;
b.l = (ymin+ymax)/2;
r.l = sqrt( sqr(a.l-xmin) + sqr(b.l-ymin) );

and now the model solves very easily.

Level values can also be set implicitly as a result of assigning bounds, since GAMS will project variable levels onto their bounds as part of executing a solve statement. For example, when a variable is bounded away from zero by a statement like Y.LO = 1; and Y is at its default level of zero, the SOLVE statement will set the level Y.L to 1.

Note: another way to formulate the model would be to minimize \(r^2\) instead of \(r\). This allows SNOPT to solve the problem even with the default starting point.

Advanced basis

GAMS automatically passes on level values and basis information from one solve to the next. Thus, when we have two solve statements in a row, with just a few changes in between SNOPT will typically need very few iterations to find an optimal solution in the second solve. For instance, when we add a second solve to the fawley.gms model from the model library:

 Model exxon /all/;
 ...
 Solve exxon maximizing profit using lp;
 Solve exxon maximizing profit using lp;

we observe the following iteration counts:

               S O L V E      S U M M A R Y

     MODEL   exxon               OBJECTIVE  profit
     TYPE    LP                  DIRECTION  MAXIMIZE
     SOLVER  SNOPT               FROM LINE  278

**** SOLVER STATUS     1 NORMAL COMPLETION
**** MODEL STATUS      1 OPTIMAL
**** OBJECTIVE VALUE             2899.2528

 RESOURCE USAGE, LIMIT          0.016      1000.000
 ITERATION COUNT, LIMIT        24         10000

.....

               S O L V E      S U M M A R Y

     MODEL   exxon               OBJECTIVE  profit
     TYPE    LP                  DIRECTION  MAXIMIZE
     SOLVER  SNOPT               FROM LINE  279

**** SOLVER STATUS     1 NORMAL COMPLETION
**** MODEL STATUS      1 OPTIMAL
**** OBJECTIVE VALUE             2899.2528

 RESOURCE USAGE, LIMIT          0.000      1000.000
 ITERATION COUNT, LIMIT         0         10000

The first solve takes 24 iterations, while the second solve needs exactly zero iterations.

Basis information is passed on using the marginals of the variables and equations. In general the rule is:

• X.M = 0: basic
• X.M \(\ne\) 0: nonbasic if level value is at bound, superbasic otherwise

A marginal value of EPS means that the numerical value of the marginal is zero, but that the status is nonbasic or superbasic. The user can specify a basis by assigning zero or nonzero values to the .M values. It is further noted that if too many .M values are zero, the basis is rejected. This happens for instance when two subsequent models are too different. This decision is made based on the value of the GAMS bratio option.

GAMS Options

The usual GAMS options (e.g. reslim, sysout) can be used to control GAMS/SNOPT. For more details, see section Controlling a Solver via GAMS Options. We highlight some of the details of this usage below for cases of special interest.

iterlim

Sets the minor iteration limit. SNOPT will stop as soon as the number of minor iterations exceeds the iteration limit, in which case the current solution will be reported.

domlim

Sets the domain error limit. Domain errors are evaluation errors in the nonlinear functions. An example of a domain error is trying to evaluate \(\sqrt{x}\) for \(x<0\). Other examples include taking logs of negative numbers, and evaluating the real power \(x^y\) for \(x < \varepsilon\) ( \(x^y\) is evaluated as \(\exp(y \log x)\)). When such an error occurs the count of domain errors is incremented: SNOPT will stop if this count exceeds the limit. If the limit has not been reached, reasonable estimates for the function (and derivatives, if requested) are returned and SNOPT continues. For example, in the case of \(\sqrt{x}, x<0\) a zero is passed back for the function value and a large value for the derivative. In many cases SNOPT will be able to recover from these domain errors, especially when they happen at some intermediate point. Nevertheless it is best to add appropriate bounds or linear constraints to ensure that these domain errors don't occur. For example, when an expression \(\log(x)\) is present in the model, add a statement like x.lo = 0.001;.

bratio

Ratio used in basis acceptance test. When a previous solution or solution estimate exists, GAMS automatically passes this solution to SNOPT so that it can reconstruct an advanced basis. When too many new (i.e. uninitialized with level and/or marginal values) variables or constraints enter the model, it may be better not to use existing basis information, but to instead crash a new basis. The bratio determines how quickly an existing basis is discarded. A value of 1.0 will discard any basis, while a value of 0.0 will retain any basis.

workfactor

By default, GAMS/SNOPT computes an estimate of the amount of workspace needed by SNOPT, and passes this workspace on to SNOPT for use in solving the model. This estimate is based on the model statistics: number of (nonlinear) equations, number of (nonlinear) variables, number of (nonlinear) nonzeroes, etc. In most cases this is sufficient to solve the model. In some rare cases SNOPT may need more memory, and the user can provide this by specifying a value of workfactor greater than 1. The computed memory estimate is multiplied by the workfactor to determine the amount of workspace made available to SNOPT for the solve.

SNOPT Options

The performance of GAMS/SNOPT is controlled by a number of parameters or "options." Each option has a default value that should be appropriate for most problems. For special situations it is possible to specify non-default values for some or all of the options via the SNOPT option file. While the content of an option file is solver-specific, the details of how to create an option file and instruct the solver to use it are not. This topic is covered in section The Solver Options File.

Note that the option file is not case sensitive. The tables below contain summary information about the SNOPT options, default values, and links to more detailed explanations.

Printing

Problem specification

Convergence tolerances

Derivative checking

Scaling

Other tolerances

QP subproblems

SQP method

Hessian approximation

Frequencies

LUSOL options

central difference interval (real): Not applicable: GAMS provides analytic derivatives ↵

Default: 6.0e-6

check frequency (integer): Controls frequency of linear constraint satisfaction test ↵

Every i^th minor iteration after the most recent basis factorization, a numerical test is made to see if the current solution \(x\) satisfies the general linear constraints (including linearized nonlinear constraints, if any). The constraints are of the form \(Ax - s = b\), where \(s\) is the set of slack variables. To perform the numerical test, the residual vector \(r = b - Ax + s\) is computed. If the largest component of \(r\) is judged to be too large, the current basis is refactorized and the basic variables are recomputed to satisfy the general constraints more accurately.

Check frequency 1 is useful for debugging purposes, but otherwise this option should not be needed.

Range: {1, ..., ∞}

Default: 60

cold start (no value): Ignore advanced basis and use CRASH procedure ↵

Requests that the CRASH procedure be used to choose an initial basis. This option takes precedence over the GAMS bratio option.

crash option (integer): Controls the basis crash algorithm ↵

Except on restarts, a CRASH procedure is used to select an initial basis from certain rows and columns of the constraint matrix \((\, A \ \ -I \,)\). The Crash option \(i\) determines which rows and columns of \(A\) are eligible initially, and how many times CRASH is called. Columns of \(-I\) are used to pad the basis where necessary.

If \(i \ge 1\), certain slacks on inequality rows are selected for the basis first. (If \(i \ge 2\), numerical values are used to exclude slacks that are close to a bound.) CRASH then makes several passes through the columns of \(A\), searching for a basis matrix that is essentially triangular. A column is assigned to "pivot" on a particular row if the column contains a suitably large element in a row that has not yet been assigned. (The pivot elements ultimately form the diagonals of the triangular basis.) For the remaining unassigned rows, slack variables are inserted to complete the basis.

By default, crash option 3 is used for linearly constrained problems and crash option 0 for problems with nonlinear constraints.

Default: auto: 0 or 3

value	meaning
0	Initial basis will be a slack basis. The initial basis contains only slack variables: \(B = I\).
1	One phase CRASH. CRASH is called once, looking for a triangular basis in all rows and columns of the matrix \(A\).
2	Two phase CRASH. CRASH is called twice (if there are nonlinear constraints). The first call looks for a triangular basis in linear rows, and the iteration proceeds with simplex iterations until the linear constraints are satisfied. The Jacobian is then evaluated for the first major iteration and CRASH is called again to find a triangular basis in the nonlinear rows (retaining the current basis for linear rows).
3	Three phase CRASH. CRASH is called up to three times (if there are nonlinear constraints). The first two calls treat linear equalities and linear inequalities separately. As before, the last call treats nonlinear rows before the first major iteration.

crash tolerance (real): Allow crash procedure to ignore small elements in eligible columns ↵

The Crash tolerance \(r\) allows the starting procedure CRASH to ignore certain small nonzeros in each column of \(A\). If \(a_{\max}\) is the largest element in column \(j\), other nonzeros \(a_{ij}\) in the column are ignored if \(|a_{ij}| \le a_{\max} \times r\). (To be meaningful, \(r\) should be in the range \(0 \le r < 1\).)

When \(r > 0.0\), the basis obtained by CRASH may not be strictly triangular, but it is likely to be nonsingular and almost triangular. The intention is to obtain a starting basis containing more columns of \(A\) and fewer (arbitrary) slacks. A feasible solution may be reached sooner on some problems.

For example, suppose the first \(m\) columns of \(A\) are the matrix shown under LU_factor_tolerance, i.e. a tridiagonal matrix with entries \(-1, 2, -1\). To help CRASH choose all \(m\) columns for the initial basis, we would specify Crash tolerance \(r\) for some value of \(r > 0.5\).

Range: [0, 1]

Default: 0.1

derivative linesearch (no value): Linesearch method (safeguarded cubic interpolation) with use of derivatives ↵

At each major iteration a linesearch is used to improve the merit function. A Derivative linesearch uses safeguarded cubic interpolation and requires both function and gradient values to compute estimates of the step size \(\alpha_k\).

difference interval (real): Not applicable: GAMS provides analytic derivatives ↵

This alters the interval \(h_1\) that is used to estimate gradients by forward differences in the following circumstances:

• In the initial ("cheap") phase of verifying the problem derivatives.
• For verifying the problem derivatives.
• For estimating missing derivatives.

In all cases, a derivative with respect to \(x_j\) is estimated by perturbing that component of \(x\) to the value \(x_j + h_1(1 + |x_j|)\), and then evaluating \(f_0(x)\) or \(f(x)\) at the perturbed point. The resulting gradient estimates should be accurate to \(O(h_1)\) unless the functions are badly scaled. Judicious alteration of \(h_1\) may sometimes lead to greater accuracy. This option has limited use in a GAMS environment as GAMS provides analytical gradients.

Default: 1.5e-8

elastic weight (real): Control for elastic mode ↵

The elastic_weight \(\omega\) determines the initial weight \(\gamma\) associated with problem NP( \(\gamma\)).

At any given major iteration \(k\), elastic mode is started if the QP subproblem is infeasible or if the QP dual variables are larger in magnitude than \(\omega(1+{\| g(x_k) \|}_2)\), where \(g\) is the objective gradient. In either case, the QP is re-solved in elastic mode with \(\gamma = \omega(1+{\| g(x_k) \|}_2)\).

Thereafter, \(\gamma\) is increased (subject to a maximum allowable value) at any point that is optimal for problem NP( \(\gamma\)), but not feasible for NP. After the \(r\)th such increase, \(\gamma=\omega 10^r(1+{\| g(x_{k1}) \|}_2)\), where \(x_{k1}\) is the iterate at which \(\gamma\) was first needed.

Default: 1.0e4

expand frequency (integer): Setting for anti-cycling mechanism ↵

This option is part of the EXPAND anti-cycling procedure [82] designed to make progress even on highly degenerate problems.

For linear models, the strategy is to force a positive step at every iteration, at the expense of violating the bounds on the variables by a small amount. Suppose that the minor feasibility tolerance is \(\delta\) and the expand frequency is \(k\). Over a period of \(k\) iterations, the tolerance actually used by SNOPT increases from \(0.5\delta\) to \(\delta\) (in steps of \(0.5\delta/k\)).

For nonlinear models, the same procedure is used for iterations in which there is only one superbasic variable. (Cycling can occur only when the current solution is at a vertex of the feasible region.) Thus, zero steps are allowed if there is more than one superbasic variable, but otherwise positive steps are enforced.

At least every k iterations, a resetting procedure eliminates any infeasible nonbasic variables. Increasing k helps to reduce the number of these slightly infeasible nonbasic variables. However, it also diminishes the freedom to choose a large pivot element (see pivot tolerance).

Range: {1, ..., ∞}

Default: 10000

factorization frequency (integer): Number of iterations between basis factorizations ↵

At most \(k\) basis changes will occur between factorizations of the basis matrix.

• With linear programs, the basis factors are usually updated every iteration. The default \(k\) is reasonable for typical problems. Smaller values (say \(k=75\) or \(k=50\)) may be more efficient on problems that are rather dense or poorly scaled.
• When the problem is nonlinear, fewer basis updates will occur as an optimum is approached. The number of iterations between basis factorizations will therefore increase. During these iterations a test is made regularly (according to the check frequency) to ensure that the general constraints are satisfied. If necessary the basis will be refactorized before the limit of \(k\) updates is reached.

By default, the frequency is set to 100 for linear models and 50 otherwise.

Range: {1, ..., ∞}

Default: auto: 100 or 50

feasible point (no value): Ignore objective function and find a feasible point ↵

The keyword feasible point means "Ignore the objective function" while finding a feasible point for the linear and nonlinear constraints. It can be used to check that the nonlinear constraints are feasible.

Default: turned off.

function precision (real): Relative accuracy with which the nonlinear functions are evaluated ↵

The relative function precision \(\epsilon_R\) is intended to be a measure of the relative accuracy with which the nonlinear functions can be computed. For example, if \(f(x)\) is computed as 1000.56789 for some relevant \(x\) and if the first 6 significant digits are known to be correct, the appropriate value for \(\epsilon_R\) would be 1.0e-6.

(Ideally the functions \(f_0(x)\) or \(f_i(x)\) should have magnitude of order 1. If all functions are substantially less than 1 in magnitude, \(\epsilon_R\) should be the absolute precision. For example, if \(f(x) =\) 1.23456789e-4 at some point and if the first 6 significant digits are known to be correct, the appropriate value for \(\epsilon_R\) would be 1.0e-10.)

• The default value of \(\epsilon_R\) is appropriate for simple analytic functions.
• In some cases the function values will be the result of extensive computation, possibly involving an iterative procedure that can provide rather few digits of precision at reasonable cost. Specifying an appropriate Function precision may lead to savings, by allowing the linesearch procedure to terminate when the difference between function values along the search direction becomes as small as the absolute error in the values.

Default: 3.00e-13

hessian frequency (integer): How often the full Hessian is reset to the identity matrix ↵

This option sets the frequency \(i\) for resetting the full memory Hessian. For example, if hessian full memory is selected and \(i\) BFGS updates have already been carried out, the Hessian approximation is reset to the identity matrix. (For certain problems, occasional resets may improve convergence, but in general they should not be necessary.)

Hessian Full Memory and Hessian Frequency = 20 have a similar effect to Hessian Limited Memory and Hessian Updates = 20, except that the latter retains the current diagonal during resets.

Default: 999999

hessian full memory (no value): Approximate Hessian is treated as a dense matrix ↵

This option selects the full storage method for storing and updating the approximate Hessian. (SNOPT uses a quasi-Newton approximation to the Hessian of the Lagrangian. A BFGS update is applied after each major iteration.)

If Hessian Full Memory is specified, the approximate Hessian is treated as a dense matrix and the BFGS updates are applied explicitly. This option is most efficient when the number of nonlinear variables \(n_1\) is not too large. In this case, the storage requirement is fixed and one can expect Q-superlinear convergence to the solution.

By default, this storage method is chosen when the number of nonlinear variables \(n_1 \le 75\).

hessian limited memory (no value): Limited-memory procedure is used to update a diagonal Hessian approximation ↵

This option selects the limited memory storage method for storing and updating the approximate Hessian. (SNOPT uses a quasi-Newton approximation to the Hessian of the Lagrangian. A BFGS update is applied after each major iteration.)

Hessian Limited Memory should be used on problems where the number of nonlinear variables \(n_1\) is large. In this case a limited-memory procedure is used to update a diagonal Hessian approximation $H_r$ a limited number of times. (Updates are accumulated as a list of vector pairs. They are discarded at regular intervals after \(H_r\) has been reset to their diagonal.)

By default, this storage method is chosen when the number of nonlinear variables \(n_1 > 75\).

hessian updates (integer): How often the limited memory Hessian is reset ↵

If hessian limited memory is selected and \(i\) BFGS updates have already been carried out, all but the diagonal elements of the accumulated updates are discarded and the updating process starts again.

Broadly speaking, the more updates stored, the better the quality of the approximate Hessian. However, the cost of each QP iteration also increases with the number of updates. The default value is likely to give a robust algorithm without significant expense, but faster convergence can sometimes be obtained with significantly fewer updates (e.g., \(i=5\)).

Default: 10

infinite bound (real): Bounds larger than this number are considered Infinity ↵

If \(r > 0\), \(r\) defines the "infinite" bound infBnd in the definition of the problem constraints. Any upper bound greater than or equal to infBnd will be regarded as plus infinity (and similarly for a lower bound less than or equal to -infBnd). If \(r \le 0\), the default value is used.

Default: 1.0e20

iterations limit (integer): Minor iteration limit ↵

The maximum number of minor iterations allowed (i.e., iterations of the simplex method or the QP algorithm), summed over all major iterations. This option, if set, overrides the GAMS iterlim specification.

Default: GAMS iterlim

linesearch tolerance (real): Controls accuracy of steplength selected ↵

This controls the accuracy with which a steplength will be located along the direction of search at each iteration. At the start of each linesearch a target directional derivative for the merit function is identified. The linesearch tolerance \(t\) determines the accuracy to which this target value is approximated.

• Larger values like \(t=0.9\) request just moderate accuracy in the linesearch.
• If the nonlinear functions are cheap to evaluate, as is usually the case for GAMS models, a more accurate search may be appropriate; try \(t=0.1\), \(0.01\) or \(0.001\). The number of major iterations might decrease.
• If the nonlinear functions are expensive to evaluate, a less accurate search may be appropriate. In the case of running under GAMS where all gradients are known, try \(t=0.99\). The number of major iterations might increase, but the total number of function evaluations may decrease enough to compensate.

Range: [0, 1]

Default: 0.9

LU complete pivoting (no value): LUSOL pivoting strategy ↵

See LU_partial_pivoting.

LU density tolerance (real): Controls when to move to a dense factorization ↵

The density tolerance \(r_1\) is used during LUSOL’s basis factorization \(B=LU\). Columns of \(L\) and rows of \(U\) are formed one at a time, and the remaining rows and columns of the basis are altered appropriately. At any stage, if the density of the remaining matrix exceeds \(r_1\), the Markowitz strategy for choosing pivots is terminated and the remaining matrix is factored by a dense \(LU\) procedure. Raising the density tolerance towards 1.0 may give slightly sparser \(LU\) factors, with a slight increase in factorization time.

See also LU_singularity_tolerance.

Range: [0, 1]

Default: 0.6

LU factor tolerance (real): Trade-off between stability and sparsity in basis factorization ↵

LU factor tolerance \(r_1\)

LU update tolerance \(r_2\)

These tolerances affect the stability and sparsity of the basis factorization \(B = LU\) during refactorization and updating, respectively. They must satisfy \(r_1\), \(r_2 \ge 1.0\). The matrix \(L\) is a product of matrices of the form

\[ \left( \begin{matrix} 1 \\ \mu & \ 1 \end{matrix} \right), \]

where the multipliers \(\mu\) satisfy \(|\mu| \le r_i\). Smaller values of \(r_i\) favor stability, while larger values favor sparsity.

• For large and relatively dense problems, smaller values of \(r_1\) (e.g. \(r_1 = 3.0\)) may give a useful improvement in stability without impairing sparsity to a serious degree.
• For certain very regular structures (e.g., band matrices) it may be necessary to reduce \(r_1\) and/or \(r_2\) in order to achieve stability. For example, if the columns of \(A\) include a submatrix of the form

  \[ \left( \begin{matrix} 2 & -1 \\ -1 & 2 & -1 \\ & -1 & 2 & -1 \\ & & \ddots& \ddots& \ddots\\ & & & -1 & 2 & -1 \\ & & & & -1 & 2 \end{matrix} \right), \]

  both \(r_1\) and \(r_2\) should be in the range \(1.0 \le r_i < 2.0\).

For linear models, the defaults are \(r_1 = 100\) and \(r_2 = 10\), while for nonlinear models both tolerances default to \(3.99\).

See also LU_update_tolerance.

Range: [1, ∞]

Default: auto

LU partial pivoting (no value): LUSOL pivoting strategy ↵

The LUSOL factorization implements a Markowitz-type search for pivots that locally minimizes fill-in subject to a threshold pivoting stability criterion. The rook pivoting and complete pivoting options are more expensive than partial pivoting but are more stable and better at revealing rank, as long as the LU_factor_tolerance is not too large (say \(r_1 < 2.0\)).

When numerical difficulties are encountered, SNOPT automatically reduces the \(LU\) tolerances toward 1.0 and switches (if necessary) to rook or complete pivoting before reverting to the default or specified options at the next refactorization. (With sysout on and system information enabled, relevant messages are output to the listing file.)

Default: yes

LU rook pivoting (no value): LUSOL pivoting strategy ↵

See LU_partial_pivoting.

LU singularity tolerance (real): Protection against ill-conditioned basis matrices ↵

The singularity tolerance \(r_2\) helps guard against ill-conditioned basis matrices. After \(B\) is refactorized, the diagonal elements of \(U\) are tested as follows: if \(|U_{jj}| \le r_2\) or \(|U_{jj}| < r_2 \max_i |U_{ij}|\), the \(j\)th column of the basis is replaced by the corresponding slack variable. (This is most likely to occur after a restart.)

See also LU_density_tolerance.

Default: 3.2e-11

LU update tolerance (real): Trade-off between stability and sparsity in basis factorization ↵

See LU_factor_tolerance for details.

Range: [1, ∞]

Default: auto

major feasibility tolerance (real): Feasibility tolerance for nonlinear constraints ↵

This tolerance \(\epsilon_r\) specifies how accurately the nonlinear constraints should be satisfied. The default value of 1.0e-6 is appropriate when the constraints are expected to have at least that accuracy.

Let rowerr be the maximum nonlinear constraint violation, normalized by the size of the solution. It is required to satisfy

\begin{equation} \tag{8} \hbox{rowerr} = \max_i\ \hbox{viol}_i / \Vert x \Vert \le\ \epsilon_r, \end{equation}

where \(\hbox{viol}_i\) is the violation of the \(i\)th nonlinear constraint.

In the GAMS/SNOPT iteration log, rowerr appears as the quantity labeled "Feasibl". If some of the problem functions are known to be of low accuracy, a larger Major feasibility tolerance may be appropriate.

Default: 1.0e-6

major iterations limit (integer): Max number of major iterations ↵

This is the maximum number of major iterations allowed. It is intended to guard against an excessive number of linearizations of the constraints. By default it is set to max(1000,m).

Default: GAMS iterlim

major optimality tolerance (real): Specifies the final accuracy of the dual variables ↵

This tolerance \(\epsilon_d\) specifies the final accuracy of the dual variables. On successful termination, SNOPT will have computed a solution \((x,s,\pi)\) such that

\begin{equation} \tag{9} \hbox{maxComp} = \max_j\ \hbox{Comp}_j / \Vert{\pi}\Vert \ \le\ \epsilon_d, \end{equation}

where \({Comp}_j\) is an estimate of the complementarity slackness for variable \(j\). The values \({Comp}_j\) are computed from the final QP solution using the reduced gradients \(d_j = g_j-\pi^T a_j\) (where \(g_j\) is the \(j\)th component of the objective gradient, \(a_j\) is the associated column of the constraint matrix \((\, A \ \ -I \,)\) and \(\pi\) is the set of QP dual variables):

\[ \hbox{Comp}_j = \left\{ \begin{array}{lcl} \phantom- d_j\min\{ x_j - l_j,1 \} & & \hbox{if $d_j \ge 0$}; \\[3pt]- d_j\min\{ u_j - x_j,1 \} & & \hbox{if $d_j < 0$.} \end{array} \right. \]

In the GAMS/SNOPT iteration log, \(maxComp\) appears as the quantity labeled "Optimal".

Default: 1.0e-6

major print level (integer): Amount of information printed during optimization (listing file) ↵

This controls the amount of output to the GAMS listing file at each major iteration. This output is only visible if the sysout option is turned on. Major print level 1 gives normal output for linear and nonlinear problems, and Major print level 11 gives additional details of the Jacobian factorization that commences each major iteration. In general, the value specified may be thought of as a binary number of the form

   Major print level   JFDXbs  

where each letter stands for a digit that is either 0 or 1 as follows:

• s single line that gives a summary of each major iteration. (This entry in JFDXbs is not strictly binary since the summary line is printed whenever \(\hbox{JFDXbs} \ge 1\)).
• b BASIS statistics, i.e., information relating to the basis matrix whenever it is refactorized. (This output is always provided if \(\hbox{JFDXbs} \ge 10\)).
• X \(x_k\), the nonlinear variables involved in the objective function or the constraints.
• D \(\pi_k\), the dual variables for the nonlinear constraints.
• F \(F(x_k)\), the values of the nonlinear constraint functions.
• J \(J(x_k)\), the Jacobian.

To obtain output of any items JFDXbs, set the corresponding digit to 1, otherwise to 0.

If J=1, the Jacobian will be output column-wise at the start of each major iteration. Column \(j\) will be preceded by the value of the corresponding variable \(x_j\) and a key to indicate whether the variable is basic, superbasic or nonbasic. (Hence if J=1, there is no reason to specify X=1 unless the objective contains more nonlinear variables than the Jacobian.) A typical line of output is

   3  1.250000D+01 BS      1  1.00000E+00      4  2.00000E+00 

which would mean that \(x_3\) is basic at value \(12.5\), and the third column of the Jacobian has elements of \(1.0\) and \(2.0\) in rows \(1\) and \(4\).

Major print level 0 suppresses most output, except for error messages.

Default: 1

major step limit (real): Limits the change in x during a linesearch ↵

This parameter \(r\) limits the change in \(x\) during a linesearch. It applies to all nonlinear problems, once a "feasible solution" or "feasible subproblem" has been found.

1. A linesearch determines a step \(\alpha\) over the range \(0 < \alpha \le \beta\), where \(\beta\) is 1 if there are nonlinear constraints, or the step to the nearest upper or lower bound on \(x\) if all the constraints are linear. Normally, the first steplength tried is \(\alpha_1 = \min(1,\beta)\).
2. In some cases, such as \(f(x) = ae^{bx}\) or \(f(x) = ax^b\), even a moderate change in the components of \(x\) can lead to floating-point overflow. The parameter \(r\) is therefore used to define a limit \( \skew{2.8}\bar\beta = r(1 + \Vert{x}\Vert) / \Vert{p}\Vert\) (where \(p\) is the search direction), and the first evaluation of \(f(x)\) is at the potentially smaller steplength \(\alpha_1 = \min(1, \skew{2.8}\bar\beta, \beta)\).
3. Wherever possible, upper and lower bounds on \(x\) should be used to prevent evaluation of nonlinear functions at meaningless points. The Major step limit provides an additional safeguard. The default value \(r = 2.0\) should not affect progress on well behaved problems, but setting \(r = 0.1\) or \(0.01\) may be helpful when rapidly varying functions are present. A "good" starting point may be required. An important application is to the class of nonlinear least-squares problems.
4. In cases where several local optima exist, specifying a small value for \(r\) may help locate an optimum near the starting point.

Default: 2.0

minor feasibility tolerance (real): Feasibility tolerance for all variables and linear constraints ↵

SNOPT tries to ensure that all variables eventually satisfy their upper and lower bounds to within this tolerance \(t\). This includes slack variables, so general linear constraints should also be satisfied to within \(t\).

Feasibility with respect to nonlinear constraints is judged by the major feasibility tolerance.

• If the bounds and linear constraints cannot be satisfied to within \(t\), the problem is declared infeasible. Let sInf be the corresponding sum of infeasibilities. If sInf is quite small, it may be appropriate to raise \(t\) by a factor of 10 or 100. Otherwise, some error in the data should be suspected.
• Nonlinear functions will be evaluated only at points that satisfy the bounds and linear constraints. If there are regions where a function is undefined, every attempt should be made to eliminate these regions from the problem. For example, if \(f(x) = \sqrt{x_1} + \log x_2\), it is essential to place lower bounds on both variables. If \(t = 1.0e-6\), the bounds \(x_1 \ge 10^{-5}\) and \(x_2 \ge 10^{-4}\) might be appropriate. (The \(\log\) singularity is more serious. In general, keep \(x\) as far away from singularities as possible.)
• If the model is scaled (see scale option), feasibility is defined in terms of the scaled problem.
• In reality, SNOPT uses \(t\) as a feasibility tolerance for satisfying the bounds on \(x\) and \(s\) in each QP subproblem. If the sum of infeasibilities cannot be reduced to zero, the QP subproblem is declared infeasible. SNOPT is then in elastic mode thereafter (with only the linearized nonlinear constraints defined to be elastic). See elastic weight for details.

Default: 1.0e-6

minor iterations limit (integer): Max number of minor iterations between linearizations of nonlinear constraints ↵

Minor iterations limit \(k\). If the number of minor iterations for the optimality phase of the QP subproblem exceeds \(k\), then all nonbasic QP variables that have not yet moved are frozen at their current values and the reduced QP is solved to optimality. Note that more than \(k\) minor iterations may be necessary to solve the reduced QP to optimality. These extra iterations are necessary to ensure that the terminated point gives a suitable direction for the linesearch. In the major iteration log, a t at the end of a line indicates that the corresponding QP was artificially terminated using the limit \(k\). Note that iterations limit defines an independent absolute limit on the total number of minor iterations (summed over all QP subproblems).

Default: 500

minor print level (integer): Amount of information printed during optimization (listing file) ↵

This controls the amount of output to the GAMS listing file during solution of the QP subproblems. It is only useful if the GAMS sysout option is turned on. The value of \(k\) has the following effect:

• \(\phantom{{}\ge}0\) No minor iteration output except error messages.
• \(\ge 1\) A single line of output each minor iteration (controlled by print frequency).
• \(\ge 10\) Basis factorization statistics generated during the periodic refactorization of the basis (see factorization frequency). Statistics for the first factorization of each major iteration are controlled by the major print level.

Default: 1

new superbasics limit (integer): Limit on new superbasics when a QP subproblem is solved ↵

This option causes early termination of the QP subproblems if the number of free variables has increased significantly since the first feasible point. If the number of new superbasics is greater than new superbasics limit the nonbasic variables that have not yet moved are frozen and the resulting smaller QP is solved to optimality. In the major iteration log, a "T" at the end of a line indicates that the QP was terminated early in this way.

Default: 99

nonderivative linesearch (no value): Linesearch method (safeguarded quadratic interpolation) without use of derivatives ↵

A nonderivative linesearch can be slightly less robust on difficult problems, and it is recommended that the default derivative linesearch be used if the functions and derivatives can be computed at approximately the same cost. If the gradients are very expensive relative to the functions, a nonderivative linesearch may give a significant decrease in computation time. In a GAMS environment derivative linesearch (the default) is more appropriate.

partial price (integer): Number of segments in partial pricing strategy ↵

This parameter sets the number of segments \(k\) using in partial pricing and is recommended for large problems that have significantly more variables than constraints. It reduces the work required for each "pricing" operation (i.e. when a nonbasic variable is selected to become superbasic).

• When \(k=1\), all columns of the constraint matrix \({(\, A \ \ -I \,)}\) are searched.
• Otherwise, \(A\) and \(I\) are partitioned to give \(k\) roughly equal segments \(A_j\), \(I_j\) ( \(j = 1\) to \(k\)). If the previous pricing search was successful on \(A_{j}\), \(I_{j}\), the next search begins on the segments \(A_{j+1}\), \(I_{j+1}\). (All subscripts here are modulo \(k\).)
• If a reduced gradient is found that is larger than some dynamic tolerance, the variable with the largest such reduced gradient (of appropriate sign) is selected to become superbasic. If nothing is found, the search continues on the next segments \(A_{j+2}\), \(I_{j+2}\), and so on.
• Partial price \(t\) (or \(t/2\) or \(t/3\)) may be appropriate for time-stage models having \(t\) time periods.

The default is 10 for linear models and 1 for nonlinear models.

Range: {1, ..., ∞}

Default: auto

penalty parameter (real): Initial penalty parameter ↵

After a QP subproblem has been solved, new estimates of the NLP solution are computed using a linesearch on the augmented Lagrangian merit function. This functions contains penalty parameters, which may be increased to ensure descent.

Default: 0

pivot tolerance (real): Used to keep the basis non-singular ↵

During solution of QP subproblems, the pivot tolerance \(r\) is used to prevent columns entering the basis if they would cause the basis to become almost singular.

• When \(x\) changes to \(x + \alpha p\) for some search direction \(p\), a ratio test is used to determine which component of \(x\) reaches an upper or lower bound first. The corresponding element of \(p\) is called the pivot element.
• Elements of \(p\) are ignored (and therefore cannot be pivot elements) if they are smaller than the pivot tolerance \(r\).
• It is common for two or more variables to reach a bound at essentially the same time. In such cases, the minor feasibility tolerance (say \(t\)) provides some freedom to maximize the pivot element and thereby improve numerical stability. Excessively small values of \(t\) should therefore not be specified.
• To a lesser extent, the expand frequency (say \(f\)) also provides some freedom to maximize the pivot element. Excessively large values of \(f\) should therefore not be specified.

Default: 3.67e-11

print frequency (integer): Number of iterations between each log line (listing file) ↵

Synonym: log_frequency

When sysout is turned on and minor print level is positive, a line of the QP iteration log will be printed on the listing file every \(k\)th minor iteration.

Range: {1, ..., ∞}

Default: 100

proximal point method (integer): Controls promimal point method used for solving linear constraints ↵

Once the linear constraints are satisfied, the proximal point method chooses a linear feasible point that is closest to \(x_0\), the initial point for the nonlinear variables. The idea is to both satisfy the linear constraints and stay close to the starting values provided for the nonlinear variables. This option is used to disable the proximal point method or to select the norm used.

Default: 1

value	meaning
0	disable PPM Do not use the proximal point method.

| 1 | one-norm
Minimize the one-norm \(||x - x_0||_1\).| | 2 | two-norm
Minimize the two-norm \(\frac{1}{2}||x - x_0||^2_2\).|

qpsolver (string): Controls method used for QP subproblems ↵

This specifies the method used to solve system (5) for the search directions in phase 2 of the QP subproblem.

• The Cholesky QP solver is the most robust, but may require a significant amount of computation if the number of superbasics is large.
• The quasi-Newton QP solver does not require the computation of the R at the start of each QP subproblem. It may be appropriate when the number of superbasics is large but relatively few major iterations are needed to reach a solution (e.g., if SNOPT is called with a warm start).
• The conjugate-gradient QP solver is appropriate for problems with many degrees of freedom (say, more than 2000 superbasics).

Default: Cholesky

value	meaning
Cholesky	full Cholesky factor QPSolver Cholesky holds the full Cholesky factor \(R\) of the reduced Hessian \(Z^THZ\). As the minor iterations proceed, the dimension of \(R\) changes with the number of superbasic variables. If the number of superbasic variables needs to increase beyond the value of reduced hessian dimension, the reduced Hessian cannot be stored and the solver switches to QPSolver CG. The Cholesky solver is reactivated if the number of superbasics stabilizes at a value less than the reduced hessian dimension.
QN	quasi-Newton method QPSolver QN solves the QP using a quasi-Newton method similar to that of MINOS. In this case, \(R\) is the factor of a quasi-Newton approximate Hessian.
CG	conjugate-gradient method QPSolver CG uses an active-set method similar to QPSolver QN, but uses the conjugate-gradient method to solve all systems involving the reduced Hessian.

reduced hessian dimension (integer): Size of Hessian matrix ↵

Synonym: hessian_dimension

This specifies that an \(i \times i\) triangular matrix \(R\) is to be available for use by the QPSolver Cholesky option (to define the reduced Hessian according to \(R^TR = Z^THZ\)). The value of \(i\) affects when QPSolver CG is activated. The default is computed internally by SNOPT, currently as \(\min\{2000, n_1 + 1\}\).

Range: {1, ..., ∞}

Default: auto

scale option (integer): Controls problem scaling ↵

Three scale options are available. By default, option 2 is used for linear models and option 1 for nonlinear models. See also scale tolerance and scale print.

Default: auto

value	meaning
0	No scaling This is recommended if it is known that \(x\) and the constraint matrix and Jacobian never have very large elements (say, larger than 100).
1	Scale linear variables Linear constraints and variables are scaled by an iterative procedure that attempts to make the matrix coefficients as close as possible to 1.0 (see Fourer [69]). This will sometimes improve the performance of the solution procedures.
2	Scale linear + nonlinear variables All constraints and variables are scaled by the iterative procedure. Also, an additional scaling is performed that takes into account columns of \((\, A \ \ -I \,)\) that are fixed or have positive lower bounds or negative upper bounds. If nonlinear constraints are present, the scales depend on the Jacobian at the first point that satisfies the linear constraints. Scale option 2 should therefore be used only if (a) a good starting point is provided, and (b) the problem is not highly nonlinear.

scale print (no value): Print scaling factors (listing file) ↵

Scale print causes the row-scales \(r(i)\) and column-scales \(c(j)\) to be printed. The scaled matrix coefficients are \({\skew3\bar a}_{ij} = a_{ij} c(j) / r(i)\), and the scaled bounds on the variables and slacks are \( {\skew2\bar l}_j = l_j / c(j),\ {\skew3\bar u}_j = u_j / c(j)\), where \(c(j) \equiv r(j-n)\) if \(j>n\).

The listing file will only show these values if the sysout option is turned on. See also scale option and scale tolerance.

scale tolerance (real): Scale tolerance ↵

The scale tolerance \(t\) affects how many passes might be needed through the constraint matrix. On each pass, the scaling procedure computes the ratio of the largest and smallest nonzero coefficients in each column:

\[\rho_j = \max_{i} |a_{ij}| / \min_i |a_{ij}| \qquad (a_{ij}\ne 0). \]

If \(\max_{j} \rho_j\) is less than \(r\) times its previous value, another scaling pass is performed to adjust the row and column scales. Raising \(r\) from 0.9 to 0.99 (say) usually increases the number of scaling passes through \(A\). At most 10 passes are made.

See also scale option and scale print.

Range: [0, 1]

Default: 0.9

solution (string): Prints SNOPT solution (listing file) ↵

This option causes the SNOPT solution file to be printed to the GAMS listing file, provided the sysout option is also turned on.

Default: turned off.

Default: NO

value	meaning
NO	Turn off printing of solution
YES	Turn on printing of solution

summary frequency (integer): Number of iterations between each log line (log file) ↵

If minor print level is positive, a line of the QP iteration log is output every \(k\)th minor iteration.

Range: {1, ..., ∞}

Default: 100

superbasics limit (integer): Maximum number of superbasics ↵

This places a limit on the storage allocated for superbasic variables. Ideally, \(i\) should be set slightly larger than the number of degrees of freedom expected at an optimal solution.

For linear programs, an optimum is normally a basic solution with no degrees of freedom. The default value of \(i\) is therefore 1. For nonlinear problems, the number of degrees of freedom is often called the "number of independent variables".

Normally, \(i\) need not be greater than \(n_1+1\), where \(n_1\) is the number of nonlinear variables. For many problems, \(i\) may be considerably smaller than \(n_1\). This will save storage if \(n_1\) is very large.

This parameter also sets the reduced hessian dimension, unless the latter is specified explicitly (and conversely). If neither parameter is specified, GAMS chooses values for both, based on problem characteristics.

Range: {1, ..., ∞}

Default: 1

suppress parameters (no value): Suppress printing of parameters (listing file) ↵

Normally SNOPT prints the option file as it is being read, and then prints a complete list of the available keywords and their final values. The suppress parameters option tells SNOPT not to print the full list. Used in conjunction with the sysout option.

system information (string): Provides additional information on the progress of the iterations (listing file) ↵

The Yes option provides additional information on the progress of the iterations, including basis repair details when ill-conditioned bases are encountered and the \(LU\) factorization parameters are strengthened.

Default: NO

value	meaning
NO	Turn off additional printing of information on progress of algorithm
YES	Turn on additional printing of information on progress of algorithm

timing level (integer): Amount of timing information (listing file) ↵

Amount of timing information written to the listing file. Used in conjunction with the sysout option.

Range: {0, ..., 3}

Default: 3

unbounded objective value (real): Determines when a problem is called unbounded ↵

This parameter sets the value \(f_{\max}\) intended to detect unboundedness in nonlinear problems. See unbounded step size for the setting of \(\alpha_{\max}\).

During a line search, \(f_0\) is evaluated at points of the form \(x + \alpha p\), where \(x\) and \(p\) are fixed and \(\alpha\) varies. If \(|f_0|\) exceeds \(f_{\max}\) or \(\alpha\) exceeds \(\alpha_{\max}\), iterations are terminated with the exit message Problem is unbounded (or badly scaled).

Unboundedness in \(x\) is best avoided by placing finite upper and lower bounds on the variables.

Default: 1.0e15

unbounded step size (real): Determines when a problem is called unbounded ↵

This parameter sets the value \(\alpha_{\max}\) used in the unboundedness test: see unbounded objective value for details.

Default: 1.0e18

verify level (integer): Finite-difference checks on the derivatives ↵

This option refers to finite-difference checks on the derivatives computed by the user-provided routines. Derivatives are checked at the first point that satisfies all bounds and linear constraints.

This option has limited use in a GAMS environment.

Default: -1

value	meaning
0	Cheap test
1	Check individual gradients
2	Check individual columns of the Jacobian
3	Combines verify level 1 and 2
-1	Derivative checking is disabled

violation limit (integer): Limit on maximum constraint violation after the linesearch ↵

This parameter \(\tau\) is used to define an absolute limit on the magnitude of the maximum constraint violation after the line search. On completion of the line search, the new iterate \(x_{k+1}\) satisfies the condition

\[ v_i(x_{k+1}) \le \tau\max\{ 1, v_i(x_0) \}, \]

where \(x_0\) is the point at which the nonlinear constraints are first evaluated and \(v_i(x)\) is the \(i\)th nonlinear constraint violation \(v_i(x) = \max( 0, l_i - F_i(x), F_i(x) - u_i )\).

The effect of this violation limit is to restrict the iterates to lie in an expanded feasible region whose size depends on the magnitude of \(\tau\). This makes it possible to keep the iterates within a region where the objective is expected to be well-defined and bounded below. If the objective is bounded below for all values of the variables, then \(\tau\) may be any large positive value.

Default: 10

warm start (no value): Use advanced basis provided by GAMS ↵

Use an advanced basis provided by GAMS. This option takes precedence over the GAMS bratio option.

The SNOPT log

When GAMS/SNOPT solves a linearly constrained problem the following log is visible on the screen:

 ==============================
 S N O P T  7.7.6    (Jan 2021)
 ==============================
 SNMEMA EXIT 100 -- finished successfully
 SNMEMA INFO 104 -- memory requirements estimated

 Nonlinear constraints       0     Linear constraints       5
 Nonlinear variables        11     Linear variables         0
 Jacobian  variables         0     Objective variables     11
 Total constraints           5     Total variables         11
 The user has defined      11   out of      11   first  derivatives

 Major Minors     Step   nObj Feasible  Optimal      Objective     nS
     0      4               1           4.3E-02 -4.7361869E+01      4   r
     1      3  1.0E+00      2           2.1E-02 -4.7526522E+01      4   r
     2      2  1.0E+00      3           3.4E-02 -4.7586643E+01      5 s
     3      2  1.0E+00      4           1.1E-02 -4.7662724E+01      6
     4      2  1.0E+00      5           4.8E-03 -4.7678618E+01      5
     5      1  1.0E+00      6           1.3E-02 -4.7701147E+01      5
     6      2  1.0E+00      7           5.2E-03 -4.7704100E+01      6
     7      2  1.0E+00      8           3.5E-04 -4.7705775E+01      5
     8      2  1.0E+00      9           5.1E-03 -4.7706132E+01      6
     9      1  1.0E+00     10           6.1E-05 -4.7706500E+01      6

 Major Minors     Step   nObj Feasible  Optimal      Objective     nS
    10      1  1.0E+00     11           6.1E-05 -4.7706513E+01      6
    11      1  1.0E+00     12           5.0E-05 -4.7706515E+01      6
    12      1  1.0E+00     13           1.6E-05 -4.7706515E+01      6
    13      1  1.0E+00     14           2.0E-06 -4.7706515E+01      6
    14      1  1.0E+00     15          (7.1E-07)-4.7706515E+01      6
 SNOPTA EXIT   0 -- finished successfully
 SNOPTA INFO   1 -- optimality conditions satisfied

 Problem name                 mixer
 No. of iterations                  30   Objective           -4.7706514831E+01
 No. of major iterations            14   Linear    obj. term  0.0000000000E+00
                                         Nonlinear obj. term -4.7706514831E+01
 User function calls (total)        15
 No. of superbasics                  6   No. of basic nonlinears             4
 No. of degenerate steps             0   Percentage                       0.00
 Max x                      11 1.6E+00   Max pi                      3 1.5E+01
 Max Primal infeas           0 0.0E+00   Max Dual infeas             4 6.1E-05

 Solution not printed
 Time for MPS input                             0.00 seconds
 Time for solving problem                       0.00 seconds
 Time for solution output                       0.00 seconds
 Time for constraint functions                  0.00 seconds
 Time for objective function                    0.00 seconds

For a nonlinearly constrained problem, the log is somewhat different:

 ==============================
 S N O P T  7.7.6    (Jan 2021)
 ==============================
 SNMEMA EXIT 100 -- finished successfully
 SNMEMA INFO 104 -- memory requirements estimated

 Nonlinear constraints      23     Linear constraints      16
 Nonlinear variables        35     Linear variables         8
 Jacobian  variables        35     Objective variables      0
 Total constraints          39     Total variables         43
 The user has defined      56   out of      56   first  derivatives

 Major Minors     Step   nCon Feasible  Optimal  MeritFunction     nS Penalty
     0     22               1  1.9E+00  1.0E+00  6.0000000E+02      3           r i
     1      3  5.0E-03      2  1.9E+00  1.0E+00 -5.8415017E+04      3 7.6E-01   rli
     2      1  5.0E-03      4  1.9E+00  1.0E+00 -2.1089442E+05      3 2.8E+00   rli
     3     19  1.0E+00      6  7.9E-01  1.7E+00 -3.1110693E+04      3 2.8E+00 sM
     4      6  1.0E+00      8  4.1E-01  4.4E+00 -3.1423080E+04      4 3.2E+00  m
     5      4  1.0E+00      9  2.1E-01  3.9E+01 -1.3134062E+04      5 3.2E+00
     6      1  1.0E+00     11  1.3E-03  2.3E+02  5.1483376E+02      5 3.2E+00  m
     7      1  1.0E+00     13  1.4E-05  2.1E+00  5.1732760E+02      5 3.2E+00  m
     8      2  1.0E+00     15  1.4E-03  3.5E+01  5.3714297E+02      4 3.2E+00  m
     9      1  1.0E+00     16  1.9E-03  4.0E+01  5.5832995E+02      4 3.2E+00

 Major Minors     Step   nCon Feasible  Optimal  MeritFunction     nS Penalty
    10      4  1.0E+00     17  1.5E-04  1.4E+01  5.6904902E+02      5 3.2E+00
    11      1  3.2E-01     20  4.9E-03  7.4E+01  5.9984627E+02      5 3.2E+00
    12      1  5.7E-01     22  7.1E-03  6.5E+01  6.3671070E+02      5 3.2E+00
    13      2  1.0E+00     23  2.6E-03  1.3E+02  6.9523972E+02      4 3.2E+00
    14      2  3.6E-01     26  5.8E-03  1.5E+02  7.3043319E+02      5 3.2E+00
[...]
    85      2  1.0E+00    281  2.6E-04  6.2E+00  1.0588496E+03      1 3.2E+00
    86      2  1.0E+00    282  2.1E-04  5.1E+00  1.0588957E+03      2 3.2E+00
    87      1  1.0E+00    283  4.9E-05  1.2E+00  1.0589180E+03      2 3.2E+00
    88      1  1.0E+00    284  5.4E-06  1.3E-01  1.0589196E+03      2 3.2E+00
    89      1  1.0E+00    285 (5.0E-07) 1.2E-02  1.0589199E+03      2 3.2E+00

 Major Minors     Step   nCon Feasible  Optimal  MeritFunction     nS Penalty
    90      1  1.0E+00    286 (7.5E-10)(1.9E-07) 1.0589199E+03      2 3.2E+00
 SNOPTA EXIT   0 -- finished successfully
 SNOPTA INFO   1 -- optimality conditions satisfied

 Problem name                 chenrad
 No. of iterations                 171   Objective            1.0589198561E+03
 No. of major iterations            90   Linear    obj. term  1.0589198561E+03
 Penalty parameter           3.188E+00   Nonlinear obj. term  0.0000000000E+00
 User function calls (total)       286
 No. of superbasics                  2   No. of basic nonlinears            29
 No. of degenerate steps             9   Percentage                       5.26
 Max x                       4 6.4E+02   Max pi                     11 2.6E+03
 Max Primal infeas          48 4.8E-07   Max Dual infeas            21 9.9E-04
 Nonlinear constraint violn    4.8E-07

 Solution not printed
 Time for MPS input                             0.00 seconds
 Time for solving problem                       0.02 seconds
 Time for solution output                       0.00 seconds
 Time for constraint functions                  0.01 seconds
 Time for objective function                    0.00 seconds

The SNOPT log shows the following columns:

Major

The current major iteration number.

Minors

The number of iterations required by both the feasibility and optimality phases of the QP subproblem. Generally, Minor will be 1 in the later iterations, since theoretical analysis predicts that the correct active set will be identified near the solution (see Section Description of the method)

Step

The step length \(\alpha\) taken along the current search direction \(p\). The variables \(x\) have just been changed to \(x+ \alpha p\). On reasonably well-behaved problems, the unit step will be taken as the solution is approached.

nObj

The number of times the nonlinear objective function has been evaluated. nObj is printed as a guide to the amount of work required for the linesearch.

nCon

The number of times SNOPT evaluated the nonlinear constraint functions.

MeritFunction

The value of the augmented Lagrangian merit function (6). This function will decrease at each iteration unless it was necessary to increase the penalty parameters (see Section Description of the method). As the solution is approached, Merit will converge to the value of the objective at the solution.

In elastic mode, the merit function is a composite function involving the constraint violations weighted by the elastic weight.

If the constraints are linear, this item is labeled Objective, the value of the objective function. It will decrease monotonically to its optimal value.

Feasible

The value of rowerr, the maximum component of the scaled nonlinear constraint residual. The solution is regarded as acceptably feasible if Feasible is less than the Major feasibility tolerance.

If the constraints are linear, all iterates are feasible and this entry is not printed.

Optimal

The value of maxgap, the maximum complementarity gap. It is an estimate of the degree of nonoptimality of the reduced costs. The solution is considered to be optimal if Optimal is less than the major optimality tolerance.

nS

The current number of superbasic variables.

Penalty

The Euclidean norm of the vector of penalty parameters used in the augmented Lagrangian merit function (not printed if the constraints are linear).

The summary line may include additional code characters that indicate what happened during the course of the iteration.

c Central differences have been used to compute the unknown components of the objective and constraint gradients. This should not happen in a GAMS environment.

d During the linesearch it was necessary to decrease the step in order to obtain a maximum constraint violation conforming to the value of violation limit.

l The norm-wise change in the variables was limited by the value of the major step limit. If this output occurs repeatedly during later iterations, it may be worthwhile increasing the value of major step limit.

i If SNOPT is not in elastic mode, an "i" signifies that the QP subproblem is infeasible. This event triggers the start of nonlinear elastic mode, which remains in effect for all subsequent iterations. Once in elastic mode, the QP subproblems are associated with the elastic problem NP( \(\gamma\)).

If SNOPT is already in elastic mode, an "i" indicates that the minimizer of the elastic subproblem does not satisfy the linearized constraints. (In this case, a feasible point for the usual QP subproblem may or may not exist.)

M An extra evaluation of the problem functions was needed to define an acceptable positive-definite quasi-Newton update to the Lagrangian Hessian. This modification is only done when there are nonlinear constraints.

m This is the same as "M" except that it was also necessary to modify the update to include an augmented Lagrangian term.

n No positive-definite BFGS update could be found. The approximate Hessian is unchanged from the previous iteration.

R The approximate Hessian has been reset by discarding all but the diagonal elements. This reset will be forced periodically by the hessian frequency and hessian updates keywords. However, it may also be necessary to reset an ill-conditioned Hessian from time to time.

r The approximate Hessian was reset after ten consecutive major iterations in which no BFGS update could be made. The diagonals of the approximate Hessian are retained if at least one update has been done since the last reset. Otherwise, the approximate Hessian is reset to the identity matrix.

s A self-scaled BFGS update was performed. This update is always used when the Hessian approximation is diagonal, and hence always follows a Hessian reset.

S This is the same as a "s" except that it was necessary to modify the self-scaled update to maintain positive definiteness.

t The minor iterations were terminated at the minor iterations limit.

T The minor iterations were terminated at the new superbasics limit.

u The QP subproblem was unbounded.

w A weak solution of the QP subproblem was found.

z The superbasics limit was reached.

Finally SNOPT prints an exit message. See Section EXIT conditions.

EXIT conditions

When the solution procedure terminates, an SNOPTA EXIT and SNOPTA INFO message is printed to summarize the final result. Here we describe each message and suggest possible courses of action.

1 -- optimality conditions satisfied

The final point seems to be a solution. This means that \(x\) is feasible (it satisfies the constraints to the accuracy requested), the reduced gradient is negligible, the reduced costs are optimal, and \(R\) is nonsingular. In all cases, some caution should be exercised. For example, if the objective value is much better than expected, SNOPT may have obtained an optimal solution to the wrong problem! Almost any item of data could have that effect if it has the wrong value. Verifying that the problem has been defined correctly is one of the more difficult tasks for a model builder.

If nonlinearities exist, one must always ask the question: could there be more than one local optimum? When the constraints are linear and the objective is known to be convex (e.g., a sum of squares) then all will be well if we are minimizing the objective: a local minimum is a global minimum in the sense that no other point has a lower function value. (However, many points could have the same objective value, particularly if the objective is largely linear.) Conversely, if we are maximizing a convex function, a local maximum cannot be expected to be global, unless there are sufficient constraints to confine the feasible region.

Similar statements could be made about nonlinear constraints defining convex or concave regions. However, the functions of a problem are more likely to be neither convex nor concave. Our advice is always to specify a starting point that is as good an estimate as possible, and to include reasonable upper and lower bounds on all variables, in order to confine the solution to the specific region of interest. We expect modelers to know something about their problem, and to make use of that knowledge as they themselves know best.

One other caution about the "Optimal solution" message. Some of the variables or slacks may lie outside their bounds more than desired, especially if scaling was requested. If sysout is on, the listing file will contain several indicators of potential issues. Max Primal infeas indicates the largest bound infeasibility and which variable is involved. If it is too large, consider restarting with a smaller minor feasibility tolerance (say 10 times smaller) and perhaps a scale option of 0.

Similarly, Max Dual infeas indicates which variable is most likely to be at a non-optimal value. Broadly speaking, if

\[ \hbox{Max Dual infeas}/\hbox{Norm of pi} = 10^{-d}, \]

then the objective function would probably change in the \(d\)th significant digit if optimization could be continued. If \(d\) seems too large, consider restarting with smaller values of major feasibility tolerance and minor feasibility tolerance.

Finally, Nonlinear constraint violn shows the maximum infeasibility for nonlinear rows. If it seems too large, consider restarting with a smaller major feasibility tolerance.

2 -- feasible point found

Occurs only if feasible point is enabled.

3 -- requested accuracy could not be achieved

If the requested accuracy could not be achieved, a feasible solution has been found, but the requested accuracy in the dual infeasibilities could not be achieved. An abnormal termination has occurred, but SNOPT is within \(10^{-2}\) of satisfying the major optimality tolerance. Check that the major optimality tolerance is not too small.

11 -- infeasible linear constraints

12 -- infeasible linear equality constraints

When the constraints are linear, the output messages are based on a relatively reliable indicator of infeasibility. Feasibility is measured with respect to the upper and lower bounds on the variables and slacks. Among all the points satisfying the general constraints \(Ax-s =0\), there is apparently no point that satisfies the bounds on \(x\) and \(s\). Violations as small as the minor feasibility tolerance are ignored, but at least one component of \(x\) or \(s\) violates a bound by more than the tolerance.

13 -- nonlinear infeasibilities minimized

14 -- linear infeasibilities minimized

15 -- infeasible linear constraints in QP subproblem

16 -- infeasible nonelastic constraints

When nonlinear constraints are present, infeasibility is much harder to recognize correctly. Even if a feasible solution exists, the current linearization of the constraints may not contain a feasible point. In an attempt to deal with this situation, when solving each QP subproblem, SNOPT is prepared to relax the bounds on the slacks associated with nonlinear rows.

If a QP subproblem proves to be infeasible or unbounded (or if the Lagrange multiplier estimates for the nonlinear constraints become large), SNOPT enters so-called "nonlinear elastic" mode. The subproblem includes the original QP objective and the sum of the infeasibilities—suitably weighted using the elastic weight parameter. In elastic mode, the nonlinear rows are made "elastic"—i.e., they are allowed to violate their specified bounds. Variables subject to elastic bounds are known as elastic variables. An elastic variable is free to violate one or both of its original upper or lower bounds. If the original problem has a feasible solution and the elastic weight is sufficiently large, a feasible point eventually will be obtained for the perturbed constraints, and optimization can continue on the subproblem. If the nonlinear problem has no feasible solution, SNOPT will tend to determine a "good" infeasible point if the elastic weight is sufficiently large. (If the elastic weight were infinite, SNOPT would locally minimize the nonlinear constraint violations subject to the linear constraints and bounds.)

Unfortunately, even though SNOPT locally minimizes the nonlinear constraint violations, there may still exist other regions in which the nonlinear constraints are satisfied. Wherever possible, nonlinear constraints should be defined in such a way that feasible points are known to exist when the constraints are linearized.

21 -- unbounded objective

22 -- constraint violation limit reached

For linear problems, unboundedness is detected by the simplex method when a nonbasic variable can apparently be increased or decreased by an arbitrary amount without causing a basic variable to violate a bound. Adding a bound on the objective will allow SNOPT to find a solution, and inspection of this solution will show the variables that can become too large due to missing restrictions.

Very rarely, the scaling of the problem could be so poor that numerical error will give an erroneous indication of unboundedness. Consider using the scale option.

For nonlinear problems, SNOPT monitors both the size of the current objective function and the size of the change in the variables at each step. If either of these is very large (see unbounded step size and unbounded objective value), the problem is terminated and declared UNBOUNDED. To avoid large function values, it may be necessary to impose bounds on some of the variables in order to keep them away from singularities in the nonlinear functions.

The second message indicates an abnormal termination while enforcing the limit on the constraint violations. This exit implies that the objective is not bounded below in the feasible region defined by expanding the bounds by the value of the violation limit.

31 -- iteration limit reached

An iteration limit was reached. Most often this is cured by increasing the GAMS iterlim option. If an SNOPT option file was used, also the iterations limit may have been set too small.

Check the iteration log to be sure that progress was being made. If so, repeat the run with higher limits. If not, consider specifying new initial values for some of the nonlinear variables.

32 -- major iteration limit reached

This indicates SNOPT was running out the limit on major iterations. This can be changed using the major iterations limit.

33 -- the superbasics limit is too small

The problem appears to be more nonlinear than anticipated. The current set of basic and superbasic variables have been optimized as much as possible and a PRICE operation is necessary to continue, but there are already as many superbasics as allowed (and no room for any more).

When increasing the superbasics limit, be aware that this also increases the reduced hessian dimension unless both options are set explicitly. This may increase the amount of memory required by SNOPT dramatically. Consider also increasing the amount memory available to SNOPT via the workfactor option, or setting a more moderate value for the reduced hessian dimension option and possibly getting slower convergence.

34 -- time limit reached

A time limit was hit. Increase the GAMS reslim option.

41 -- current point cannot be improved

The algorithm could not find a better solution although optimality was not achieved within the optimality tolerance. Possibly scaling can lead to better function values and derivatives. Raising the major optimality tolerance will probably make this message go away. Try better scaling, better bounds or a better starting point.

42 -- singular basis

The first factorization attempt found the basis to be structurally or numerically singular. (Some diagonals of the triangular matrix \(U\) were deemed too small.) The associated variables were replaced by slacks and the modified basis refactorized, but singularity persisted. Try better scaling, better bounds or a better starting point.

43 -- cannot satisfy the general constraints

The basic variables \(x_B\) have been recomputed, given the present values of the superbasic and nonbasic variables. A step of "iterative refinement" has also been applied to increase the accuracy of \(x_B\), but a row check has revealed that the resulting solution does not satisfy the QP constraints \(Ax - s = b\) sufficiently well. Try better scaling, better bounds or a better starting point.

44 -- ill-conditioned null-space basis

During computation of the reduced Hessian \(Z^THZ\), some column(s) of \(Z\) continued to contain very large values. Try better scaling, better bounds or a better starting point.

51 -- incorrect objective derivatives

52 -- incorrect constraint derivatives

The derivatives are not deemed to be correct. This message should not occur using a GAMS model without external equations.

53 -- irregular or badly scaled problem functions

The problem function can be badly scaled in case of DNLPs (which SNOPT was not designed for) or because of evaluations at singular points.

61 -- undefined function at the first feasible point

62 -- undefined function at the initial point

SNOPT was unable to proceed because the functions are undefined at the initial point or the first feasible point. Try to add better bounds or linear equations such that non-linear functions can be evaluated or use a better starting point.

63 -- unable to proceed into undefined region

Repeated attempts to move into a region where the functions are not defined resulted in the change in variables being unacceptably small. At the final point, it appears that the only way to decrease the merit function is to move into a region where the problem functions are not defined.

Try to add better bounds or linear equations such that non-linear functions can be evaluated or use a better starting point.

71 -- terminated during function evaluation

72 -- terminated from monitor routine

These messages indicate trouble evaluating the non-linear functions or derivatives. Usually these errors show a "Function evaluation error limit" message.

81 -- work arrays must have at least 500 elements

82 -- not enough character storage

83 -- not enough integer storage

84 -- not enough real storage

Increase the memory allocation by using workfactor > 1.