NAG FL Interface
e04srf (handle_​solve_​ssqp)

Note: this routine uses optional parameters to define choices in the problem specification and in the details of the algorithm. If you wish to use default settings for all of the optional parameters, you need only read Sections 1 to 10 of this document. If, however, you wish to reset some or all of the settings please refer to Section 11 for a detailed description of the algorithm and to Section 12 for a detailed description of the specification of the optional parameters.
Note: please be advised that this routine is classed as ‘experimental’ and its interface may be developed further in the future. Please see Section 4 in How to Use the NAG Library for further information.
Settings help

FL Name Style:


FL Specification Language:


1 Purpose

e04srf is a solver from the NAG optimization modelling suite for large-scale Nonlinear Programming (NLP) problems based on an active-set Sequential Quadratic Programming (SQP) method, using limited-memory quasi-Newton approximations to the Hessian of the Lagrangian.

2 Specification

Fortran Interface
Subroutine e04srf ( handle, objfun, objgrd, confun, congrd, hess, monit, nvar, x, nnzu, u, rinfo, stats, iuser, ruser, cpuser, ifail)
Integer, Intent (In) :: nvar, nnzu
Integer, Intent (Inout) :: iuser(*), ifail
Real (Kind=nag_wp), Intent (Inout) :: x(nvar), u(nnzu), ruser(*)
Real (Kind=nag_wp), Intent (Out) :: rinfo(100), stats(100)
Type (c_ptr), Intent (In) :: handle, cpuser
External :: objfun, objgrd, confun, congrd, hess, monit
C Header Interface
#include <nag.h>
void  e04srf_ (void **handle,
void (NAG_CALL *objfun)(const Integer *nvar, const double x[], double *fx, Integer *inform, Integer iuser[], double ruser[], void **cpuser),
void (NAG_CALL *objgrd)(const Integer *nvar, const double x[], const Integer *nnzfd, double fdx[], Integer *inform, Integer iuser[], double ruser[], void **cpuser),
void (NAG_CALL *confun)(const Integer *nvar, const double x[], const Integer *ncnln, double gx[], Integer *inform, Integer iuser[], double ruser[], void **cpuser),
void (NAG_CALL *congrd)(const Integer *nvar, const double x[], const Integer *nnzgd, double gdx[], Integer *inform, Integer iuser[], double ruser[], void **cpuser),
void (NAG_CALL *hess)(const Integer *nvar, const double x[], const Integer *ncnln, const Integer *idf, const double *sigma, const double lambda[], const Integer *nnzh, double hx[], Integer *inform, Integer iuser[], double ruser[], void **cpuser),
void (NAG_CALL *monit)(const Integer *nvar, const double x[], const Integer *nnzu, const double u[], Integer *inform, const double rinfo[], const double stats[], Integer iuser[], double ruser[], void **cpuser),
const Integer *nvar, double x[], const Integer *nnzu, double u[], double rinfo[], double stats[], Integer iuser[], double ruser[], void **cpuser, Integer *ifail)
The routine may be called by the names e04srf or nagf_opt_handle_solve_ssqp.

3 Description

e04srf is designed to minimize (maximize) a linear or nonlinear function possibly subject to a set of constraints including bounds on the variables and sparse linear and nonlinear constraints. It is suitable for large-scale general nonlinear programming (NLP) problems, as well as for linearly constrained optimization. It is particularly 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, there is no limit on the number of degrees of freedom.
e04srf solves a nonlinear programming problem in the following form
minimize xn f(x) subject to lgg(x)ug, 12 xTQix + riTx + si0 ,   1imQ, lBBxuB, lxxux, (1)
where
e04srf serves as a solver for problems stored as a handle. The handle points to an internal data structure which defines the problem and serves as a means of communication for routines in the NAG optimization modelling suite. First, the problem handle is initialized, typically by calling e04raf. Then some of the components of the model need to be defined by calling one or more of the following routines based on the component type:
Once the problem is fully described, the handle may be passed to the solver e04srf. When the handle is no longer needed, e04rzf should be called to destroy it and deallocate the memory held within. There are other ways that the problem model can be defined and/or updated, see Section 3.1 in the E04 Chapter Introduction for more details about the NAG optimization modelling suite.
The algorithm behaviour can be modified by various optional parameters (see Section 12) which can be set by e04zmf and e04zpf anytime between the initialization of the handle and a call to the solver. Once the solver has finished, options or the model may be modified for the next solve. The solver may be called repeatedly with various starting points and/or optional parameters. Option getter e04znf can be called to retrieve the current value of any option.
The optional parameter Task may be used to switch the problem to maximization.
In certain situations the solver can be sped up by warm starting, see Section 9.4 and optional parameter SSQP Start Type for details.
Several options may have a significant impact on the performance of the solver. Even if the defaults were chosen to suit the majority of problems, it is recommended that you experiment in order to find the most suitable set of options for a particular problem, see Sections 11 and 12 for further details.

3.1 Note regarding nonlinear functions

In the presence of a nonlinear objective function or constraint function, the structure (sparsity pattern) of these functions are defined by e04rgf and e04rkf respectively; the function values and the first derivatives (gradients) are computed at requested points by the corresponding supplied functions objfun, objgrd, confun and congrd. Ideally, the first derivatives should be known and coded but e04srf estimates any missing ones by finite differences, see SSQP Estimate Derivatives.
The solver can strongly benefit from the correct identification of linear variables; these are the variables which appear only in linear expressions in the whole optimization model. The solver cannot autodetect if any variable present in nonlinear constraints or the nonlinear objective is linear thus there is a mechanism to declare the variable as linear by e04rcf and it is highly recommended to do so whenever possible. See Section 3.1 in e04rcf for a detailed explanation and examples. The main advantage is that linear variables have constant first partial derivatives and must be zero in second derivatives, thus the internal approximation of the Hessian doesn't need to take them into account which speeds up the solver.

3.2 Alternative solvers

e04srf is a generic solver which can solve a great variety of problems, however, a dedicated solver might work better in the following cases. See Section 4 in the E04 Chapter Introduction.

4 References

Dantzig G B (1963) Linear Programming and Extensions Princeton University Press
Eldersveld S K (1991) Large-scale sequential quadratic programming algorithms PhD Thesis Department of Operations Research, Stanford University, Stanford
Fourer R (1982) Solving staircase linear programs by the simplex method Math. Programming 23 274–313
Gill P E, Murray W and Saunders M A (2002) SNOPT: An SQP Algorithm for Large-scale Constrained Optimization 12 979–1006 SIAM J. Optim.
Gill P E, Murray W and Saunders M A (2015) Users' guide for SNOPT 7.5: Software for large-scale linear nonlinear programming Report SNOPT Department of Mathematics, University of California, San Diego https://www.ccom.ucsd.edu/~peg/papers/sndoc7.pdf
Gill P E, Murray W and Saunders M A (2016) Users' guide for SQOPT 7.5: Software for large-scale linear and quadratic programming Report SQOPT Department of Mathematics, University of California, San Diego https://www.ccom.ucsd.edu/~peg/papers/sqdoc7.pdf
Gill P E, Murray W, Saunders M A and Wright M H (1986) Users' guide for NPSOL (Version 4.0): a Fortran package for nonlinear programming Report SOL 86-2 Department of Operations Research, Stanford University
Gill P E, Murray W, Saunders M A and Wright M H (1987) Maintaining LU factors of a general sparse matrix Linear Algebra and its Applics. 88/89 239–270
Gill P E, Murray W, Saunders M A and Wright M H (1992) Some theoretical properties of an augmented Lagrangian merit function Advances in Optimization and Parallel Computing (ed P M Pardalos) 101–128 North Holland
Hock W and Schittkowski K (1981) Test Examples for Nonlinear Programming Codes. Lecture Notes in Economics and Mathematical Systems 187 Springer–Verlag
Murtagh B A and Saunders M A (1978) Large-scale linearly constrained optimization 14 41–72 Math. Programming
Murtagh B A and Saunders M A (1982) A projected Lagrangian algorithm and its implementation for sparse nonlinear constraints Math. Program. Stud. 16 84–118
Murtagh B A and Saunders M A (1995) MINOS 5.4 users' guide Report SOL 83-20R Department of Operations Research, Stanford University

5 Arguments

1: handle Type (c_ptr) Input
On entry: the handle to the problem. It needs to be initialized (e.g., by e04raf) and to hold a problem formulation compatible with e04srf. It must not be changed between calls to the NAG optimization modelling suite.
2: objfun Subroutine, supplied by the NAG Library or the user. External Procedure
objfun must calculate the value of the nonlinear objective function f(x) at a specified point x. If there is no nonlinear objective (i.e., e04rgf was not called), objfun will never be called by e04srf and objfun may be the dummy routine e04srv included in the NAG Library.
The specification of objfun is:
Fortran Interface
Subroutine objfun ( nvar, x, fx, inform, iuser, ruser, cpuser)
Integer, Intent (In) :: nvar
Integer, Intent (Inout) :: inform, iuser(*)
Real (Kind=nag_wp), Intent (In) :: x(nvar)
Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
Real (Kind=nag_wp), Intent (Out) :: fx
Type (c_ptr), Intent (In) :: cpuser
C Header Interface
void  objfun (const Integer *nvar, const double x[], double *fx, Integer *inform, Integer iuser[], double ruser[], void **cpuser)
1: nvar Integer Input
On entry: n, the current number of decision variables x in the model.
2: x(nvar) Real (Kind=nag_wp) array Input
On entry: the vector x of variable values at which the objective function is to be evaluated.
3: fx Real (Kind=nag_wp) Output
On exit: the value of the objective function at x.
4: inform Integer Input/Output
On entry: a non-negative value.
On exit: may be used to indicate that the function cannot be evaluated at the requested point x by setting inform<0. The algorithm will try to recover if the objective cannot be evaluated; if recovery is not possible it will stop with ifail=21 or ifail=25.
5: iuser(*) Integer array User Workspace
6: ruser(*) Real (Kind=nag_wp) array User Workspace
7: cpuser Type (c_ptr) User Workspace
objfun is called with the arguments iuser, ruser and cpuser as supplied to e04srf. You should use the arrays iuser and ruser, and the data handle cpuser to supply information to objfun.
objfun must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which e04srf is called. Arguments denoted as Input must not be changed by this procedure.
3: objgrd Subroutine, supplied by the NAG Library or the user. External Procedure
objgrd must calculate the values of the nonlinear objective function gradient f x at a specified point x. If there is no nonlinear objective (i.e., e04rgf was not called), objgrd will never be called by e04srf and objgrd may be the dummy subroutine e04srw (included in the NAG Library.)
If the optional parameter SSQP Estimate Derivatives=YES or ONLYOBJ, then after returning from objgrd the gradient vector is checked for missing entries which you have not supplied. Missing entries will be estimated using a finite difference method, see optional parameter SSQP Estimate Derivatives description for more details.
Note: objgrd should not return floating-point NaN (Not a Number) or infinity values, if your code inadvertently returns any NaNs or infinities, e04srf is likely to terminate with ifail=21, 24 or 25.
The specification of objgrd is:
Fortran Interface
Subroutine objgrd ( nvar, x, nnzfd, fdx, inform, iuser, ruser, cpuser)
Integer, Intent (In) :: nvar, nnzfd
Integer, Intent (Inout) :: inform, iuser(*)
Real (Kind=nag_wp), Intent (In) :: x(nvar)
Real (Kind=nag_wp), Intent (Inout) :: fdx(nnzfd), ruser(*)
Type (c_ptr), Intent (In) :: cpuser
C Header Interface
void  objgrd (const Integer *nvar, const double x[], const Integer *nnzfd, double fdx[], Integer *inform, Integer iuser[], double ruser[], void **cpuser)
1: nvar Integer Input
On entry: n, the current number of decision variables x in the model.
2: x(nvar) Real (Kind=nag_wp) array Input
On entry: the vector x of variable values at which the objective function gradient is to be evaluated.
3: nnzfd Integer Input
On entry: the number of nonzero elements in the sparse gradient vector of the objective function, as was set in a previous call to e04rgf.
4: fdx(nnzfd) Real (Kind=nag_wp) array Input/Output
On entry: the elements should only be assigned and not referenced.
On exit: the values of the nonzero elements in the sparse gradient vector of the objective function, in the order specified by idxfd in a previous call to e04rgf. fdx(i) will store the gradient element f xj , where j=idxfd(i).
5: inform Integer Input/Output
On entry: a non-negative value.
If inform=100, then this is a special presolve call. This call is used to fill-in the constant elements of the objective gradient, e.g., when there are linear variables in the model. Note that this call may occur with an iterate x that may not be feasible with respect to the linear constraints but it will always be feasible with respect to the variable bounds. Furthermore, this call only requires to return the constant elements of the gradient, f xi, where f(x) is linear in xi. If you do not provide these elements, subsequent calls to objfun with inform=0 may occur in the vicinity of x with the purpose of estimating the missing derivatives using finite differences.
On exit: may be used to inform that the gradient cannot be evaluated at the requested point x by setting inform<0. The algorithm will try to recover if the gradient cannot be evaluated; if recovery is not possible it will stop with ifail=21 or ifail=25.
6: iuser(*) Integer array User Workspace
7: ruser(*) Real (Kind=nag_wp) array User Workspace
8: cpuser Type (c_ptr) User Workspace
objgrd is called with the arguments iuser, ruser and cpuser as supplied to e04srf. You should use the arrays iuser and ruser, and the data handle cpuser to supply information to objgrd.
objgrd must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which e04srf is called. Arguments denoted as Input must not be changed by this procedure.
4: confun Subroutine, supplied by the NAG Library or the user. External Procedure
confun must calculate the values of the nonlinear constraints g(x) at a specified point x. If there are no nonlinear constraints then confun will never be called by e04srf and it may be the dummy subroutine e04srx included in the NAG Library.
Note: confun should not return floating-point NaN (Not a Number) or infinity values, if your code inadvertently returns any NaNs or infinities, e04srf is likely to terminate with ifail=21, 24 or 25.
The specification of confun is:
Fortran Interface
Subroutine confun ( nvar, x, ncnln, gx, inform, iuser, ruser, cpuser)
Integer, Intent (In) :: nvar, ncnln
Integer, Intent (Inout) :: inform, iuser(*)
Real (Kind=nag_wp), Intent (In) :: x(nvar)
Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
Real (Kind=nag_wp), Intent (Out) :: gx(ncnln)
Type (c_ptr), Intent (In) :: cpuser
C Header Interface
void  confun (const Integer *nvar, const double x[], const Integer *ncnln, double gx[], Integer *inform, Integer iuser[], double ruser[], void **cpuser)
1: nvar Integer Input
On entry: n, the current number of decision variables x in the model.
2: x(nvar) Real (Kind=nag_wp) array Input
On entry: the vector x of variable values at which the constraint functions are to be evaluated.
3: ncnln Integer Input
On entry: mg, the number of nonlinear constraints, as specified in an earlier call to e04rkf.
4: gx(ncnln) Real (Kind=nag_wp) array Output
On exit: the mg values of the nonlinear constraint functions at x.
5: inform Integer Input/Output
On entry: a non-negative value.
If inform=100, then this is a special presolve call. This call is used to quickly identify if the problem is infeasible, e.g., due to fixed variables. Note that this call can occur with an iterate x that may not be feasible with respect to the linear constraints but it will always be feasible with respect to the variable bounds.
On exit: may be used to inform that the constraint values cannot be evaluated at the requested point x by setting inform<0. The algorithm will try to recover if the gradient cannot be evaluated; if recovery is not possible it will stop with ifail=21 or ifail=25.
6: iuser(*) Integer array User Workspace
7: ruser(*) Real (Kind=nag_wp) array User Workspace
8: cpuser Type (c_ptr) User Workspace
confun is called with the arguments iuser, ruser and cpuser as supplied to e04srf. You should use the arrays iuser and ruser, and the data handle cpuser to supply information to confun.
confun must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which e04srf is called. Arguments denoted as Input must not be changed by this procedure.
5: congrd Subroutine, supplied by the NAG Library or the user. External Procedure
congrd must calculate the Jacobian of the nonlinear constraints g x at a specified point x. If there are no nonlinear constraints, congrd will never be called by e04srf and congrd may be the dummy subroutine e04sry included in the NAG Library.
If the optional parameter SSQP Estimate Derivatives=YES or ONLYCON, then after returning from congrd the Jacobian is checked for missing entries which you have not supplied. Missing entries will be estimated using a finite difference method, see optional parameter SSQP Estimate Derivatives description for more details.
Note: congrd should not return floating-point NaN (Not a Number) or infinity values, if your code inadvertently returns any NaNs or infinities, e04srf is likely to terminate with ifail=21, 24 or 25.
The specification of congrd is:
Fortran Interface
Subroutine congrd ( nvar, x, nnzgd, gdx, inform, iuser, ruser, cpuser)
Integer, Intent (In) :: nvar, nnzgd
Integer, Intent (Inout) :: inform, iuser(*)
Real (Kind=nag_wp), Intent (In) :: x(nvar)
Real (Kind=nag_wp), Intent (Inout) :: gdx(nnzgd), ruser(*)
Type (c_ptr), Intent (In) :: cpuser
C Header Interface
void  congrd (const Integer *nvar, const double x[], const Integer *nnzgd, double gdx[], Integer *inform, Integer iuser[], double ruser[], void **cpuser)
1: nvar Integer Input
On entry: n, the current number of decision variables x in the model.
2: x(nvar) Real (Kind=nag_wp) array Input
On entry: the vector x of variable values at which the Jacobian of the constraint functions is to be evaluated.
3: nnzgd Integer Input
On entry: is the number of nonzero elements in the sparse Jacobian of the constraint functions, as was set in a previous call to e04rkf.
4: gdx(nnzgd) Real (Kind=nag_wp) array Input/Output
On entry: the elements should only be assigned and not referenced.
On exit: the nonzero values of the Jacobian of the nonlinear constraints, in the order specified by irowgd and icolgd in an earlier call to e04rkf. gdx(i) will be the gradient gj xk , where j=irowgd(i) and k=icolgd(i).
5: inform Integer Input/Output
On entry: a non-negative value.
If inform=100, then this is a special presolve call. This call is used to fill-in the constant elements of the Jacobian, e.g., when there are linear variables in the model. Note that this call may occur with an iterate x that may not be feasible with respect to the linear constraints but it will always be feasible with respect to the variable bounds. Furthermore, this call only requires to return the constant elements of the Jacobian, gj xi , where gj(x) is linear in xi. If the user does not provide these elements, subsequent calls to confun with inform=0 may occur in the vicinity of x with the purpose of estimating the missing derivatives using finite differences.
On exit: may be used to inform that the Jacobian of the constraints cannot be evaluated at the requested point x by setting inform<0. The algorithm will try to recover if the gradient cannot be evaluated; if recovery is not possible it will stop with ifail=21 or ifail=25.
6: iuser(*) Integer array User Workspace
7: ruser(*) Real (Kind=nag_wp) array User Workspace
8: cpuser Type (c_ptr) User Workspace
congrd is called with the arguments iuser, ruser and cpuser as supplied to e04srf. You should use the arrays iuser and ruser, and the data handle cpuser to supply information to congrd.
congrd must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which e04srf is called. Arguments denoted as Input must not be changed by this procedure.
6: hess Subroutine, supplied by the NAG Library or the user. External Procedure
hess is reserved for future releases of the NAG Library which will allow calculation of the second derivatives. It will never be called in the current implementation and hess may be the dummy subroutine e04srz included in the NAG Library.
The specification of hess is:
Fortran Interface
Subroutine hess ( nvar, x, ncnln, idf, sigma, lambda, nnzh, hx, inform, iuser, ruser, cpuser)
Integer, Intent (In) :: nvar, ncnln, idf, nnzh
Integer, Intent (Inout) :: inform, iuser(*)
Real (Kind=nag_wp), Intent (In) :: x(nvar), sigma, lambda(ncnln)
Real (Kind=nag_wp), Intent (Inout) :: hx(nnzh), ruser(*)
Type (c_ptr), Intent (In) :: cpuser
C Header Interface
void  hess (const Integer *nvar, const double x[], const Integer *ncnln, const Integer *idf, const double *sigma, const double lambda[], const Integer *nnzh, double hx[], Integer *inform, Integer iuser[], double ruser[], void **cpuser)
1: nvar Integer Input
On entry: n, the current number of decision variables x in the model.
2: x(nvar) Real (Kind=nag_wp) array Input
On entry: the vector x of variable values at which the Hessian functions are to be evaluated.
3: ncnln Integer Input
On entry: mg, the number of nonlinear constraints, as specified in an earlier call to e04rkf.
4: idf Integer Input
On entry: specifies the quantities to be computed in hx.
idf=−1
The values of the Hessian of the Lagrangian will be computed in hx. This will be the case if e04rlf has been called with idf of the same value.
idf=0
The values of the Hessian of the objective function will be computed in hx. This will be the case if e04rlf has been called with idf of the same value.
idf>0
The values of the Hessian of the constraint function with index idf will be computed in hx. This will be the case if e04rlf has been called with idf of the same value.
5: sigma Real (Kind=nag_wp) Input
On entry: if idf=−1, the value of the σ quantity in the definition of the Hessian of the Lagrangian. Otherwise, sigma should not be referenced.
6: lambda(ncnln) Real (Kind=nag_wp) array Input
On entry: if idf=−1, the values of the λi quantities in the definition of the Hessian of the Lagrangian. Otherwise, lambda should not be referenced.
7: nnzh Integer Input
On entry: the number of nonzero elements in the Hessian to be computed.
8: hx(nnzh) Real (Kind=nag_wp) array Input/Output
On entry: the elements should only be assigned and not referenced.
On exit: the nonzero values of the requested Hessian evaluated at x. For each value of idf, the ordering of nonzeros must follow the sparsity structure registered in the handle by earlier calls to e04rlf through the arguments irowh and icolh.
9: inform Integer Input/Output
On entry: a non-negative value.
On exit: must be set to a value describing the action to be taken by the solver on return from hess. Specifically, if the value is negative the solution of the current problem will terminate immediately with ifail=21 or ifail=25; otherwise, computations will continue.
10: iuser(*) Integer array User Workspace
11: ruser(*) Real (Kind=nag_wp) array User Workspace
12: cpuser Type (c_ptr) User Workspace
hess is called with the arguments iuser, ruser and cpuser as supplied to e04srf. You should use the arrays iuser and ruser, and the data handle cpuser to supply information to hess.
hess must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which e04srf is called. Arguments denoted as Input must not be changed by this procedure.
7: monit Subroutine, supplied by the NAG Library or the user. External Procedure
monit is provided to enable you to monitor the progress of the optimization and optionally to terminate the solver early if necessary. It is invoked at the end of every ith major iteration where i is given by the SSQP Monitor Frequency (the default is 0, monit is not called).
monit may be the dummy subroutine e04sru included in the NAG Library.
The specification of monit is:
Fortran Interface
Subroutine monit ( nvar, x, nnzu, u, inform, rinfo, stats, iuser, ruser, cpuser)
Integer, Intent (In) :: nvar, nnzu
Integer, Intent (Inout) :: inform, iuser(*)
Real (Kind=nag_wp), Intent (In) :: x(nvar), u(nnzu), rinfo(100), stats(100)
Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
Type (c_ptr), Intent (In) :: cpuser
C Header Interface
void  monit (const Integer *nvar, const double x[], const Integer *nnzu, const double u[], Integer *inform, const double rinfo[], const double stats[], Integer iuser[], double ruser[], void **cpuser)
1: nvar Integer Input
On entry: n, the current number of decision variables x in the model.
2: x(nvar) Real (Kind=nag_wp) array Input
On entry: xi, the values of the decision variables x at the current iteration.
3: nnzu Integer Input
On entry: the dimension of array u.
4: u(nnzu) Real (Kind=nag_wp) array Input
On entry: if nnzu>0, u holds the values of Lagrange multipliers (dual variables, see Section 11.3) for the constraints at the current iteration. See Section 3.1 in e04stf for layout information.
5: inform Integer Input/Output
On entry: a non-negative value.
On exit: may be used to request the solver to stop immediately. Specifically, if inform<0 the solver will terminate immediately with ifail=20; otherwise, the solver will proceed normally.
6: rinfo(100) Real (Kind=nag_wp) array Input
On entry: error measures and various indicators at the end of the current iteration as described in rinfo.
7: stats(100) Real (Kind=nag_wp) array Input
On entry: solver statistics at the end of the current iteration as described in stats.
8: iuser(*) Integer array User Workspace
9: ruser(*) Real (Kind=nag_wp) array User Workspace
10: cpuser Type (c_ptr) User Workspace
monit is called with the arguments iuser, ruser and cpuser as supplied to e04srf. You should use the arrays iuser and ruser, and the data handle cpuser to supply information to monit.
monit must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which e04srf is called. Arguments denoted as Input must not be changed by this procedure.
8: nvar Integer Input
On entry: n, the current number of decision variables x in the model.
9: x(nvar) Real (Kind=nag_wp) array Input/Output
On entry: x0, the initial estimates of the variables, x.
On exit: the final values of the variables, x.
10: nnzu Integer Input
On entry: the number of Lagrange multipliers in array u.
If nnzu=0, u will not be referenced; otherwise, it needs to match the dimension q as explained in Section 3.1 in e04stf.
Constraints:
  • nnzu0;
  • if nnzu>0, nnzu=q.
11: u(nnzu) Real (Kind=nag_wp) array Input/Output
Note: if nnzu>0, u holds Lagrange multipliers (dual variables) for the constraints. See Section 3.1 in e04stf for layout information. If nnzu=0, u will not be referenced.
On entry: optionally provides the initial estimates of Lagrange multipliers, (z0,π0). These values are only referenced if a warm start is requested (Section 9.4). If there are no initial estimates available, then set to zero.
On exit: the final value of Lagrange multipliers (z*,π*), see Section 11.3.
12: rinfo(100) Real (Kind=nag_wp) array Output
On exit: error measures and various indicators at the end of the final iteration as given in the table below:
1 Objective function value f(x).
2 Constraint violation (primal infeasibility), see Section 7.
3 Dual infeasibility.
47 Reserved for future use.
8 Step size for the primal variables, see Iteration log in Section 9.1.
9 Estimate of the condition number of RTR , itself an estimate of ZTHZ , the reduced Hessian of the Lagrangian. See Section 9.1.
10 Largest nonlinear constraint violation.
11 Largest nonlinear constraint violation, normalized using 1+ xk, with xk the latest iterate. See (12).
12 Sum of the infeasibilities of constraints that lie outside one of their bounds by more than the optional parameter SSQP Minor Feasibility Tol before the solution is unscaled.
13 Flag indicating whether the current iterate satisfies the first-order Kuhn–Tucker conditions for optimality:
j Description
0 none of the conditions are satisfied;
1 condition on primal infeasibility satisfied, rinfo(2) is within SSQP Major Feasibility Tol tolerance;
2 condition on dual infeasibility satisfied, rinfo(3) is within SSQP Minor Feasibility Tol tolerance; and
3 both conditions are satisfied.
14 The size of the penalty parameter D of (8). See Section 11.4.
15 Flag indicating whether the current iterate is feasible with respect to the linear constraints, if any. A nonzero value indicates feasibility.
16100 Reserved for future use.
13: stats(100) Real (Kind=nag_wp) array Output
On exit: solver statistics at the end of the final iteration as given in the table below:
1 Number of major iterations.
2 Total number of minor iterations.
3 Number of minor iterations performed on the last major iteration.
4 Number of Hessian evaluations.
5 Number of gradient evaluations.
6
Number of infeasibilities of constraints that lie outside one of their bounds by more than the optional parameter SSQP Minor Feasibility Tol before the solution is unscaled.
If any linear constraints are infeasible, it contains the number of variables and linear constraints lying outside their upper or lower bounds. The nonlinear constraints are not evaluated.
Otherwise, it contains the number of components of g(x) lying outside their bounds by more than the optional parameter SSQP Minor Feasibility Tol. Again this is before the solution is unscaled.
7 Total time spent in presolver (including user-supplied function calls).
8 Total time spent in solver (including user-supplied function calls).
9 Total time spent evaluating constraint function.
10 Total time spent evaluating constraint Jacobian.
11 Total time spent evaluating objective function.
12 Total time spent evaluating objective gradient.
13 Total time spent evaluating Hessian.
1418 Reserved for future use.
19 Number of function evaluations.
20 Number of constraint function evaluations.
21 Number of constraint Jacobian evaluations.
2224 Reserved for future use.
25 Number of superbasic variables.
26 Number of column swaps between B and S.
27 Number of elements in the LU factorization of B.
2829 Reserved for future use.
30 Number of slack variables. See Sections 11.1 and 9.2.
31100 Reserved for future use.
14: iuser(*) Integer array User Workspace
15: ruser(*) Real (Kind=nag_wp) array User Workspace
16: cpuser Type (c_ptr) User Workspace
iuser, ruser and cpuser are not used by e04srf, but are passed directly to objfun, objgrd, confun, congrd, hess and monit and may be used to pass information to these routines. If you do not need to reference cpuser, it should be initialized to c_null_ptr.
17: ifail Integer Input/Output
On entry: ifail must be set to 0, −1 or 1 to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of 0 causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of −1 means that an error message is printed while a value of 1 means that it is not.
If halting is not appropriate, the value −1 or 1 is recommended. If message printing is undesirable, then the value 1 is recommended. Otherwise, the value −1 is recommended since useful values can be provided in some output arguments even when ifail0 on exit. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry ifail=0 or −1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
Note: in some cases e04srf may return useful information.
ifail=1
The supplied handle does not define a valid handle to the data structure for the NAG optimization modelling suite. It has not been properly initialized or it has been corrupted.
ifail=2
The problem is already being solved.
This solver does not support the model defined in the handle.
ifail=4
On entry, nvar=value, expected value=value.
Constraint: nvar must match the current number of variables of the model in the handle.
ifail=5
On entry, nnzu=value.
Constraint: nnzu=value or 0.
On entry, nnzu=value.
Constraint: no constraints present, so nnzu must be 0.
ifail=7
The dummy confun routine was called but the problem requires these values. Please provide a proper confun routine.
The dummy congrd routine was called but the problem requires these derivatives. Please provide a proper congrd routine.
The dummy hess routine was called but the problem requires these derivatives. Please provide a proper hess routine.
The dummy objfun routine was called but the problem requires these values. Please provide a proper objfun routine.
The dummy objgrd routine was called but the problem requires these derivatives. Please provide a proper objgrd routine.
ifail=20
User requested termination during a monitoring step. inform was set to a negative value in monit.
ifail=21
The current starting point is unusable.
Either inform was set to a negative value within the user-supplied functions objfun, objgrd, confun, congrd or hess, or an Infinity or NaN was detected in values returned from them. Note that the user-provided starting point is made bound-constrain feasible (if applicable) before calling any of the user-supplied functions.
See also ifail=25.
ifail=22
Maximum number of iterations reached.
Either the SSQP Iteration Limit, SSQP Major Iteration Limit or SSQP Minor Iteration Limit limit was exceeded before a solution with the requested accuracy could be found. Check the iteration log to be sure that progress was being made. If so, rerun the problem after increasing these optional parameters.
ifail=23
The solver terminated after the maximum time allowed was exhausted.
Maximum number of seconds exceeded. Use optional parameter Time Limit to change the limit.
ifail=24
The solver was terminated because no further progress could be achieved.
This can indicate that numerical difficulties have been encountered and little or no progress was made. It could also happen if the problem has been solved to the best numerical accuracy possible given the current scaling. Additionally, several other circumstances could lead to this error:
  1. 1.objfun or confun could be returning accurate function values but objgrd or congrd inaccurate gradients (or vice versa). This is the most likely cause. Study the comments given for ifail=26 as well as those given for Verify Derivatives and SSQP Estimate Derivatives. Ensure that the coding is correct.
  2. 2.The function and gradient values could be consistent, but their precision could be too low. For example, accidental use of a low precision data type when a higher precision was intended would lead to a relative function precision of about 10−6 instead of something like 10−15 . The default SSQP Major Optimality Tol of 2×10−6 would need to be raised to about 10−3 for optimality to be declared (at a rather suboptimal point). Of course, it is better to revise the function coding to obtain as much precision as economically possible.
  3. 3.If function values are obtained from an expensive iterative process, they may be accurate to rather few significant figures, and gradients will probably not be available. One should specify but even then, if t is as large as 10−5 or 10−6 (only 5 or 6 significant figures), the same exit condition may occur. At present, the only remedy is to increase the accuracy of the function calculation.
  4. 4.An LU factorization of the basis has just been obtained and used to recompute the basic variables xB , given the present values of the superbasic and nonbasic variables. A step of ‘iterative refinement’ has also been applied to increase the accuracy of xB . However, a row check has revealed that the resulting solution does not satisfy the current constraints Ax-s=0 sufficiently well. This probably means that the current basis is very ill-conditioned. If there are some linear constraints and variables, try SSQP Scale Option=VARS AND LC if scaling has not yet been used.
  5. 5.The first factorization attempt will have found the basis to be structurally or numerically singular (some diagonals of the triangular matrix U were respectively zero or smaller than a certain tolerance.) The associated variables are replaced by slacks and the modified basis is refactorized, but singularity persists. This should mean that the problem is badly scaled.
ifail=25
Invalid value detected in user function and recovery failed.
Either inform was set to a negative value within the user-supplied functions objfun, objgrd, confun, congrd or hess, or an Infinity or NaN was detected in values returned from them. This will occur in the case when one or more of the problem functions are
  • undefined in a region around the current iterate and no recovery is possible;
  • undefined at the first feasible point; or
  • undefined when checking derivatives.
In these cases, it may be convenient to add a linear constraint block to exclude the undefined region. See also ifail=21.
Additionally, other circumstances could lead to this error:
  • The solver is detecting irregular or badly scaled problem functions. See optional parameter SSQP Function Precision.
  • The solver iterates are in a region where the QP Hessian is numerically indefinite and the optimization process was stopped.
ifail=26
User-provided derivatives are likely to be incorrect.
This error indicates that the user-supplied gradient vector fdx or Jacobian matrix gdx have at least one value that disagrees remarkably with its associated forward difference estimate (the relative difference between the computed and estimated values is 1.0 or more). This exit is a safeguard since e04srf will usually fail to make progress when the computed gradients are seriously inaccurate. In the process it may expend considerable effort before terminating with ifail=24.
Check the function and Jacobian computation very carefully in objgrd and congrd. A simple omission could explain everything. If a component is very large, then give serious thought to scaling the function or the nonlinear variables. If you feel certain that the computed Jacobian is correct (and that the forward-difference estimate is, therefore, wrong), you can specify Verify Derivatives=NO to prevent individual elements from being checked. However, the optimization procedure may have difficulties or even fail.
ifail=27
Maximum number of superbasic variables reached.
Increase value of the optional parameter SSQP Superbasics Limit.
This error may mean that 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 pricing operation (where a nonbasic variable is selected to become superbasic) is necessary to continue, but it cannot continue as the number of superbasic variables has already reached the limit specified by the optional parameter SSQP Superbasics Limit. In general, raise the limit by a reasonable amount, bearing in mind the storage needed for the reduced Hessian.
ifail=50
Problem was solved to an acceptable level, full accuracy was not achieved.
This indicates that the algorithm found a feasible solution, but the requested accuracy in the dual infeasibilities could not be achieved. The solution found is within 10−2 of satisfying the SSQP Major Feasibility Tol. This may happen if the SSQP Major Feasibility Tol tolerance is too small for the current problem, or if the problem is badly scaled.
ifail=51
The problem was found to be infeasible during preprocessing.
One or more of the constraints (or its part after preprocessing) violates the constraints.
ifail=53
The problem seems to be infeasible, the algorithm was stopped.
When the constraints are linear, this message is 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 (see (5) and (6) in Section 11.2), there is apparently no point that satisfies the bounds on x and s . Violations as small as the SSQP Minor Feasibility Tol are ignored, but at least one component of x or s violates a bound by more than the tolerance.
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, e04srf 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), e04srf enters a so-called ‘nonlinear elastic’ mode. The subproblem includes the original QP objective and the sum of the infeasibilities — suitably weighted. In elastic mode, some of the bounds on the nonlinear rows are ‘elastic’, i.e., they are allowed to violate their specific 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 will eventually be obtained for the perturbed constraints, and optimization can continue on the subproblem. If the nonlinear problem has no feasible solution, e04srf will tend to determine a ‘good’ infeasible point if the elastic weight is sufficiently large. If the elastic weight were infinite, e04srf would locally minimize the nonlinear constraint violations subject to the linear constraints and bounds.
Unfortunately, even though e04srf 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.
ifail=54
The problem seems to be unbounded and the algorithm was stopped.
This indicates that the solver found a direction where the objective value decreases arbitrarily. In the case of maximizing, the direction found increases arbitrarily the objective value. Cases where this happens are when a linear objective is provided with no constraints, or when the constraints defined do not bound the search direction.
For linear problems, unboundedness is detected by the simplex method when a nonbasic variable can be increased or decreased by an arbitrary amount without causing a basic variable to violate a bound. Consider adding an upper or lower bound to the variable. Also, examine the constraints that have nonzeros in the associated column, to see if they have been formulated as intended.
Very rarely, the scaling of the problem could be so poor that numerical error will give an erroneous indication of unboundedness. Consider using the optional parameter SSQP Scale Option.
For nonlinear problems, e04srf 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, 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.
ifail=90
Memory limit is too small.
e04srf has terminated because the LU factorization exhausted the allocated memory limit. This indicates that the memory allocation heuristic underestimated the requirements. Please contact NAG.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

If the value of the optional parameter SSQP Major Optimality Tol is set to 10-d and ifail=0 on exit, then the final value of f(x) should have approximately d correct significant digits.
e04srf returns with ifail=0 if the iterates have converged to a point x that satisfies the first-order Kuhn–Tucker (see Section 9.1.6) conditions to the accuracy requested by the optional parameter SSQP Major Optimality Tol, i.e., the projected gradient and active constraint residuals are negligible at x.
You should check whether the following four conditions are satisfied:
  1. (i)the final value of rgNorm (see Section 9.1.6) is significantly less than that at the starting point;
  2. (ii)during the final major iterations, the values of Step and Minors (see Section 9.1.5) are both one;
  3. (iii)the last few values of both rgNorm and SumInf (see Section 9.1.6) become small at a fast linear rate; and
  4. (iv)condHz (see Section 9.1.5) is small.
If all these conditions hold, x is almost certainly a local minimum of (1).
A cautionary note about ‘Optimal solutions’. Some of the variables or slacks may lie outside their bounds more than desired, especially if scaling was requested. Max Primal infeas (printed at the end of the solve when Print Level and Monitoring Level4 or Print Level and Monitoring Level4) refers to the largest bound infeasibility and which variable is involved. If it is too large, consider restarting with a smaller SSQP Minor Feasibility Tol (say 10 times smaller) and perhaps SSQP Scale Option=NONE.
Similarly, Max Dual infeas (also printed when Print Level or Monitoring Level4) indicates which variable is most likely to be at a nonoptimal value. Broadly speaking, if
Max Dual infeas / Max pi = 10-d ,  
where Max pi (printed at the end of a solve when Print Level and Monitoring Level1 or Print Level and Monitoring Level1) indicates the largest dual variable in magnitude, 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 a smaller SSQP Major Optimality Tol.
Finally, Nonlinear constraint violn (printed at the end of a solve when Print Level and Monitoring Level1 or Print Level and Monitoring Level1) reports the maximum infeasibility for nonlinear rows. If it seems too large, consider restarting with a smaller SSQP Major Feasibility Tol.

8 Parallelism and Performance

Background information to multithreading can be found in the Multithreading documentation.
e04srf is not threaded in any implementation.

9 Further Comments

9.1 Description of the Printed Output

The solver can print information to give an overview of the problem and of the progress of the computation. The output may be sent to two independent streams (files) which are set by optional parameters Print File and Monitoring File. Optional parameters Print Level, Monitoring Level, Print Solution and Print Options determine the exposed level of detail. This allows, for example, a detailed log file to be generated while the condensed information is displayed on the screen.
By default (Print File=6, Print Level=2), the following sections are printed to the standard output:

9.1.1 Header

The header is a message indicating the start of the solver. It should look like:
-------------------------------------------------------------------
 E04SR, SQP method for large-scale nonlinear optimization problems
-------------------------------------------------------------------

9.1.2 Optional parameters list

The list shows all options of the solver, each displayed on one line. The output contains the option name, its current value and an indicator for how it was set. The options unchanged from the default setting are noted by ‘d’, options you set are noted by ‘U’, and options reset by the solver are noted by ‘S’. Note that the output format is compatible with the file format expected by e04zpf. The output might look as follows:
Begin of Options
  Print File                    =                   6     * d
  Print Level                   =                   2     * U
  Print Options                 =                 Yes     * d
  Print Solution                =                 Yes     * U
  Monitoring File               =                   9     * U
  Monitoring Level              =                   3     * U

  Infinite Bound Size           =         1.00000E+20     * d
  Task                          =            Minimize     * d
  Stats Time                    =                  No     * d
  Time Limit                    =         1.00000E+06     * d
  Verify Derivatives            =                 Yes     * U

  ...
  Ssqp Estimate Derivatives     =                  No     * d
  Ssqp Function Precision       =         3.00000E-13     * d
  Ssqp Finite Diff Ctrl Interval=         6.69433E-05     * d
  Ssqp Finite Diff Interval     =         5.47723E-07     * d
  Ssqp Hessian                  =                Auto     * d
  Ssqp Hessian Updates          =                  10     * d
  Ssqp Iteration Limit          =               10000     * d
  Ssqp Major Feasibility Tol    =         1.00001E-06     * d
  Ssqp Major Iteration Limit    =                1000     * d
  ...
End of Options

9.1.3 Problem statistics

Problem statistics
If Print Level2, statistics on the problem are printed. It should look as follows:
Problem Statistics
 No of variables                  6
   linear                         4
   nonlinear                      2
   free (unconstrained)           2
   bounded                        2
 No of lin. constraints           0
   nonzeroes                      0
 No of quad.constraints           9
 No of nln. constraints          12
   nonzeroes                     15
 Objective function       Quadratic
Note: the list is showing that out of the six variables, four are detected as present in linear expressions (of the objective and constraint functions), while the remaining two are present nonlinearly.
It is followed by information regarding warm start (see Section 9.4) and statistics on the Jacobian matrix, the output may look similar to:
 Warm start information loaded successfully from handle.
       Handle [WARM START BASIS] data origin: solver
       Handle [SLACK] data origin           : solver

 Matrix statistics
               Total      Normal        Free       Fixed     Bounded
 Rows              4           1           1           2           0
 Columns           4           2           2           0           0

 The user has defined       5   out of       5   constraint gradients.
 The user has defined       4   out of       4   objective  gradients.

9.1.4 Verification of derivatives

If Verify Derivatives=CHEAP, then only a ‘cheap’ test will be performed on the user-supplied nonlinear objective gradient and constraint Jacobian. The output may be similar to:
 Cheap test of user-supplied problem derivatives...
 The constraint gradients seem to be OK.
 -->  The largest discrepancy was    2.52E-07  in constraint     3 (NLC:      2)
 The objective  gradients seem to be OK.
 Gradient projected in one direction   6.70906812736E-01
 Difference approximation              6.70906898365E-01
If Verify Derivatives=YES, then the user-supplied nonlinear objective gradient and constraint Jacobian are individually verified. The output might look as follows:
 Verification of user-supplied problem derivatives.
 The constraint gradients seem to be OK.
 -->  The largest discrepancy was    2.52E-07  in constraint     3 (NLC:      2)
 The objective  gradients seem to be OK.
 Gradient projected in one direction   6.70906812736E-01
 Difference approximation              6.70906898365E-01
 Column       x(j)        dx(j)    Element no.    Row        Derivative    Difference approxn
      1 -7.06220960E-02  5.86E-07         1         2   -1.41244192E-01   -1.41243605E-01  ok
      1 -7.06220960E-02  5.86E-07           Objective    2.68365409E+00    2.68365467E+00  ok
      2  1.41244914E+00  1.32E-06         2         2    2.82489828E+00    2.82489960E+00  ok
                                          3         3    1.12714152E+01    1.12714310E+01  ok
      2  1.41244914E+00  1.32E-06           Objective    2.68365409E+00    2.68365441E+00  ok
...
 -->  The largest relative error was    1.29E-06   in row     3 (NLC:      2),  column     2
      4  objective gradients out of     1  through     4  seem to be OK.
 -->  The largest relative error was    3.59E-07   in column     2

9.1.5 Major Iteration Log

This section describes the output to unit Print File or Monitoring File if Print Level​ or ​Monitoring Level>1. One line of information is output every major iteration. A heading is printed before the first such line containing the items described below.
Label Description
Itns is the cumulative number of minor iterations.
Major is the current major iteration number.
Minors is the number of iterations required by both the feasibility and optimality phases of the QP subproblem. Generally, Minors will be 1 in the later iterations, since theoretical analysis predicts that the correct active set will be identified near the solution, see Section 11.
Step is the step length α taken along the current search direction p . The variables x have just been changed to x + α p . On reasonably well-behaved problems, the unit step will be taken as the solution is approached.
nCon the number of times objfun or confun have been called to evaluate the nonlinear problem functions. Evaluations needed for the estimation of the derivatives by finite differences are not included. nCon is printed as a guide to the amount of work required for the linesearch.
Feasible is the value of vmax from (12), the maximum component of the scaled nonlinear constraint residual (see optional parameter SSQP Major Feasibility Tol). The solution is regarded as acceptably feasible if Feasible is less than the SSQP Major Feasibility Tol. In this case, the entry is contained in parentheses.
If the constraints are linear, all iterates are feasible and this entry is not printed.
Optimal is the value of cmax from (13), the maximum complementary gap (see optional parameter SSQP Major Optimality Tol). It is an estimate of the degree of nonoptimality of the reduced costs. Both Feasible and Optimal are small in the neighbourhood of a solution.
MeritFunction is the value of the augmented Lagrangian merit function from (7). This function will decrease at each iteration unless it was necessary to increase the penalty parameters, see Section 11.4. As the solution is approached, MeritFunction 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 labelled Objective, the value of the objective function. It will decrease monotonically to its optimal value.
L+U is the number of nonzeros representing the basis factors L and U on completion of the QP subproblem.
If nonlinear constraints are present, the basis factorization B = LU is computed at the start of the first minor iteration. At this stage, L+U = lenL+lenU , where lenL is the number of subdiagonal elements in the columns of a lower triangular matrix and lenU is the number of diagonal and superdiagonal elements in the rows of an upper-triangular matrix, see Section 9.1.9.
As columns of B are replaced during the minor iterations, L+U may fluctuate up or down but, in general, will tend to increase. As the solution is approached and the minor iterations decrease towards zero, L+U will reflect the number of nonzeros in the LU factors at the start of the QP subproblem.
BSwap is the number of columns of the basis matrix B that were swapped with columns of S to improve the condition of B . The swaps are determined by an LU factorization of the rectangular matrix BS = (B S) T with stability being favoured more than sparsity.
nS is the current number of superbasic variables.
condHz is an estimate of the condition number of RTR , itself an estimate of ZTHZ , the reduced Hessian of the Lagrangian. The condition number is the square of the ratio of the largest and smallest diagonals of the upper triangular matrix R , this being a lower bound on the condition number of RTR . condHz gives a rough indication of whether or not the optimization procedure is having difficulty. If ε is the relative machine precision being used, the SQP algorithm will make slow progress if condHz becomes as large as ε-1/2 108 , and will probably fail to find a better solution if condHz reaches ε-3/4 1012 .
To guard against high values of condHz, attention should be given to the scaling of the variables and the constraints. In some cases it may be necessary to add upper or lower bounds to certain variables to keep them a reasonable distance from singularities in the nonlinear functions or their derivatives.
Penalty is the Euclidean norm of the vector of penalty parameters used in the augmented Lagrangian merit function (not printed if there are no nonlinear constraints).
The summary line may include additional code characters that indicate what happened during the course of the major iteration. These will be under the ...Flags column.
Flag Description
+ indicates that the log is from a major iteration.
c central differences have been used to compute the unknown components of the objective and constraint gradients. A switch to central differences is made if either the linesearch gives a small step, or x is close to being optimal. In some cases, it may be necessary to re-solve the QP subproblem with the central difference gradient and Jacobian.
d during the linesearch it was necessary to decrease the step in order to obtain a maximum constraint violation.
D inform<0 on return from objfun, indicating that the linesearch needed to be done with a smaller value of the step length αk of equation (9).
l the norm wise change in the variables was limited by a safeguard upper bound.
i if e04srf 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 (11) described in Section 11.5.
If e04srf 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 SSQP Hessian Updates keyword. 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 used when the Hessian approximation is diagonal, and hence always follows a Hessian reset.
t the minor iterations were terminated because of the SSQP Minor Iteration Limit.
T the minor iterations were terminated because the limit on new superbasics was reached.
u the QP subproblem was unbounded.
w a weak solution of the QP subproblem was found.
z the SSQP Superbasics Limit was reached.

9.1.6 Minor Iteration Log

Is only printed if Print Level​ or ​Monitoring Level 3 , one line of information is output every k th minor iteration, where k is the specified SSQP Print Frequency, by default k=1. A heading is printed before the first such line following a basis factorization. The heading contains the items described below. In this description, a pricing operation is the process by which a nonbasic variable is selected to become superbasic (in addition to those already in the superbasic set). The selected variable is denoted by jq. Variable jq often becomes basic immediately. Otherwise, it remains superbasic unless it reaches its opposite bound and returns to the nonbasic set.
Label Description
Itn the current iteration number.
LPmult or QPmult is the reduced cost (or reduced gradient) of the variable jq selected by the pricing procedure at the start of the present iteration. Algebraically, the reduced gradient is dj = gj - πT aj for j = jq , where gj is the gradient of the current objective function, π is the vector of dual variables for the QP subproblem, and aj is the j th column of ( A -I ) .
Note that the reduced cost is the 1-norm of the reduced-gradient vector at the start of the iteration, just after the pricing procedure.
LPstep or QPstep is the step length α taken along the current search direction p . The variables x have just been changed to x + α p . Write Step to stand for LPStep or QPStep, depending on the problem. If a variable is made superbasic during the current iteration ( +SBS>0 ), Step will be the step to the nearest bound. During Phase 2, the step can be greater than 1 only if the reduced Hessian is not positive definite.
nInf is the number of infeasibilities after the present iteration. This number will not increase unless the iterations are in elastic mode.
SumInf is the sum of infeasibilities after the present iteration, if nInf>0. The value usually decreases at each nonzero Step, but if it decreases by 2 or more, SumInf may occasionally increase.
rgNorm is the norm of the reduced-gradient vector at the start of the iteration. (It is the norm of the vector with elements dj for variables j in the superbasic set.) During Phase 2 this norm will be approximately zero after a unit step. (The heading is not printed if the problem is linear.)
LPobjective or QPobjective the QP objective function after the present iteration. In elastic mode, the heading is changed to Elastic QPobj. In either case, the value printed decreases monotonically.
+SBS is the variable jq selected by the pricing operation to be added to the superbasic set.
-SBS is the superbasic variable chosen to become nonbasic.
-BS is the basis variable removed (if any) to become nonbasic.
Pivot if column aq replaces the r th column of the basis B , Pivot is the r th element of a vector y satisfying By=aq . Wherever possible, Step is chosen to avoid extremely small values of Pivot (since they cause the basis to be nearly singular).
L+U is the number of nonzeros representing the basis factors L and U . Immediately after a basis factorization B=LU , L+U is lenL+lenU, the number of subdiagonal elements in the columns of a lower triangular matrix and the number of diagonal and superdiagonal elements in the rows of an upper-triangular matrix. Further nonzeros are added to L when various columns of B are later replaced. As columns of B are replaced, the matrix U is maintained explicitly (in sparse form). The value of L will steadily increase, whereas the value of U may fluctuate up or down. Thus the value of L+U may fluctuate up or down (in general, it will tend to increase).
ncp is the number of compressions required to recover storage in the data structure for U . This includes the number of compressions needed during the previous basis factorization.
nS is the current number of superbasic variables. (The heading is not printed if the problem is linear.)
condHz see Section 9.1.5. (The heading is not printed if the problem is linear.)

9.1.7 Summary

Once the solver finishes, a detailed summary is produced:
 Problem                      Minimize
 Max x                       2 1.4E+00   Max pi                      2 1.9E+01
 Nonlinear constraint violn    5.2E-08
 
 -------------------------------------------------------------------------------
 Status: converged, an optimal solution was found
 -------------------------------------------------------------------------------
 
 Value of the objective               1.90012E+00
 Primal infeasibility                 5.24140E-08
 Dual infeasibility                   5.19087E-06
 Penalty parameter                    9.04109E+00
 Mayor iterations                              11
 Minor iterations                              11
 Objective function calls                      18
 Objective gradient calls                      10
 Constraint function calls                     16
 Constraint Jacobian calls                     10
 Hessian evaluations                            0
 -------------------------------------------------------------------------------
It starts with a status line of the overall result, followed by the final objective value as well as the primal and dual measures of infeasibility and the final value of the penalty parameter. If Print Level​ and ​Monitoring Level2, it will additionally report major and minor iteration count, objective or constraint function calls. Optionally, if Stats Time=YES or CPU, some statistics on timings is displayed. It might look as follows:
 Timing                                        Seconds         Average
   Total                                         18.32                  100.0%
   Total time spent in pre-solve                  1.09                    5.9%
   Total time spent in solver                    17.23                   94.1%
   Time in objective function                     0.31            0.02    1.7%
   Time in objective gradient                     1.24            0.12    6.8%
   Time in constraint function                    0.45            0.03    2.5%
   Time in constraint Jacobian                    2.68            0.27   14.6%
   Time in Hessian callback                       0.00            0.00    0.0%

9.1.8 Solution

If Print Solution=YES, the values of the primal variables x are printed, furthermore if the problem is constrained, also the dual variables (Lagrangian multipliers) u and the constraint bounds are reported. It might look as follows:
 Primal variables:
   idx   Lower bound       Value       Upper bound
     1       -inf      -7.06221E-02         inf
     2       -inf       1.41245E+00         inf
...
 
 Box bounds dual variables:
   idx   Lower bound       Value       Upper bound       Value
     1       -inf       0.00000E+00         inf       2.92732E-06
     2       -inf       0.00000E+00         inf       5.31212E-08
...
 
 Quadratic constraints dual variables:
   idx   Lower bound       Value       Upper bound       Value
     1       -inf       0.00000E+00    0.00000E+00   -0.00000E+00
...
 
 Nonlinear constraints dual variables:
   idx   Lower bound       Value       Upper bound       Value
     1   2.00000E+00    0.00000E+00    2.00000E+00    1.90001E+01
     2   4.00000E+00    5.00000E+00    4.00000E+00    0.00000E+00
...

9.1.9 Basis Factorization Statistics

If Print Level​ or ​Monitoring Level 4 , the first seven items listed below are printed whenever the basis B or the rectangular matrix BS = (B S) T is factorized before finding the solution of the next QP subproblem. Note that BS may be factorized at the start of just some of the major iterations. It is immediately followed by a factorization of B itself.
Gaussian elimination is used to compute a sparse LU factorization of B or BS , where PLPT and PUQ are lower and upper triangular matrices, for some permutation matrices P and Q .
If Print Level​ or ​Monitoring Level 4 , the same items are printed during the QP solution whenever the current B is factorized.
A heading is printed containing the items described below.
Label Description
Factor the number of factorizations since the start of the run.
Demand a code giving the reason for the present factorization, as follows:
Code Meaning
0 First LU factorization.
1 The number of updates reached the allowed limit.
2 The nonzeros in the updated factors have increased significantly.
7 Not enough storage to update factors.
10 Row residuals are too large.
11 Ill-conditioning has caused inconsistent results.
Itn is the current minor iteration number.
Nonlin is the number of nonlinear variables in the current basis B .
Linear is the number of linear variables in B .
Slacks is the number of slack variables in B .
B, BR, BS or BT factorize is the type of LU factorization.
B periodic factorization of the basis B .
BR more careful rank-revealing factorization of B using threshold rook pivoting. This occurs mainly at the start, if the first basis factors seem singular or ill-conditioned. Followed by a normal B factorize.
BS BS is factorized to choose a well-conditioned B from the current (B S) . Followed by a normal B factorize.
BT same as BS except the current B is tried first and accepted if it appears to be not much more ill-conditioned than after the previous BS factorize.
m is the number of rows in B or BS .
n is the number of columns in B or BS . Preceded by ‘=’ or ‘>’ respectively.
Elems is the number of nonzero elements in B or BS .
Amax is the largest nonzero in B or BS .
Density is the percentage nonzero density of B or BS .
Merit/MerRP/MerCP Merit is the average Markowitz merit count for the elements chosen to be the diagonals of PUQ . Each merit count is defined to be (c-1) (r-1) where c and r are the number of nonzeros in the column and row containing the element at the time it is selected to be the next diagonal. Merit is the average of n such quantities. It gives an indication of how much work was required to preserve sparsity during the factorization. Depending on the LU pivoting strategy, this heading may change to MerCP or MerRP.
lenL is the number of nonzeros in L .
L+U is the number of nonzeros representing the basis factors L and U . Immediately after a basis factorization B=LU , this is lenL+lenU, the number of subdiagonal elements in the columns of a lower triangular matrix and the number of diagonal and superdiagonal elements in the rows of an upper-triangular matrix. Further nonzeros are added to L when various columns of B are later replaced. As columns of B are replaced, the matrix U is maintained explicitly (in sparse form). The value of L will steadily increase, whereas the value of U may fluctuate up or down. Thus the value of L+U may fluctuate up or down (in general, it will tend to increase).
Cmpressns is the number of times the data structure holding the partially factored matrix needed to be compressed to recover unused storage. Ideally this number should be zero. If it is more than 3 or 4, the amount of workspace available to e04srf should be increased for efficiency.
Incres is the percentage increase in the number of nonzeros in L and U relative to the number of nonzeros in B or BS .
Utri is the number of triangular rows of B or BS at the top of U .
lenU the number of nonzeros in U , including its diagonals.
Ltol is the largest subdiagonal element allowed in L .
Umax the maximum nonzero element in U .
Ugrwth is the ratio Umax/Amax , which ideally should not be substantially larger than 10.0 or 100.0.
As long as Lmax is not large (say 5.0 or less), max(Amax,Umax) / DUmin gives an estimate of the condition number B . If this is extremely large, the basis is nearly singular. Slacks are used to replace suspect columns of B and the modified basis is refactored.
Ltri is the number of triangular columns of B or BS at the left of L .
dense1 is the number of columns remaining when the density of the basis matrix being factorized reached 0.3.
Lmax is the actual maximum subdiagonal element in L (bounded by Ltol).
Akmax is the largest nonzero generated at any stage of the LU factorization. (Values much larger than Amax indicate instability.)
Agrwth is the ratio Akmax/Amax . Values much larger than 100 (say) indicate instability.
bump is the size of the block to be factorized nontrivially after the triangular rows and columns of B or BS have been removed.
dense2 is the number of columns remaining when the density of the basis matrix being factorized reached 0.6. (The Markowitz pivot strategy searches fewer columns at that stage.)
DUmax is the largest diagonal element of PUQ .
DUmin is the smallest diagonal element of PUQ .
condU the ratio DUmax/DUmin , which estimates the condition number of U (and of B if Ltol is less than 5.0, say).

9.2 Accessing Slack Variables

e04srf reformulates any inequality constraint in the model into an equality constraint. This is done implicitly by adding artificial slack variables to the problem, see Section 11.1.
e04srf stores in the handle under the label ‘SLACK’ the final value of the slack variables. It also retrieves this information from the handle when a warm start is requested using the optional parameter SSQP Start Type=WARM, see Section 9.4.
Slack variables can be stored or retrieved from the handle with e04rxf using cmdstr='SLACK' and lrarr=stats(30), where stats(30) is the correct number of slack variables and can be retrieved from a previous call to e04srf.
Note: any entry of the slack variables array with an index greater than the amount of defined variables and constraints will refer to the artificial slack variable associated with the objective function, this may be the case when linear variables are identified in the objective function.

9.3 Retrieving and Storing a Basis

A basis refers to a partitioning of the primal and slack variables. This partitioning plays a fundamental role in the underlying LP and QP algorithms of e04srf, Section 11.3 provides a brief description.
e04srf stores in the handle under the label ‘BASIS’ (or ‘WARM START BASIS’) the final state of the primal and slack variables. It also retrieves this information from the handle when a warm start is requested using optional parameter SSQP Start Type=WARM, see Section 9.4.
The stored integer array is of length nvar+stats(30) and the values describe the state of the primal variables x and the slacks as follows:
BASIS(i) State of variable i Usual value Note
-1 Ignored None of the enabled variables can have this value. It indicates that the basis element is ignored by the solver. Ignored elements occur when a primal variable or constraint is marked by you as disabled (see e04tcf) or presolve dropped a fixed constraint or primal variable.
0 Nonbasic lower bound
1 Nonbasic upper bound
2 Superbasic between lower and
upper bounds
3 Basic between lower and
upper bounds
4 Nonbasic lower bound Variable is ignored during the crash procedure and its associated slack variable needs to be set to the lower bound.
5 Nonbasic upper bound Variable is ignored during the crash procedure and its associated slack variable needs to be set to the upper bound.
The basis can be stored or retrieved from the handle with e04rwf using cmdstr='BASIS' or 'WARM START BASIS' and liarr=nvar+stats(30), where stats(30) is the correct number of slack variables and can be retrieved from a previous call to e04srf.
Note: basic and superbasic variables may be outside their bounds by as much as the SSQP Minor Feasibility Tol tolerance. If scaling is specified (using optional parameter SSQP Scale Option), the feasibility tolerance applies to the variables of the scaled problem. In this case, the variables of the original problem may be as much as 0.1 outside their bounds, but this is unlikely unless the problem is very badly scaled. Check the ‘Primal infeasibility’ line printed at the end of the solve, see Section 9.1.7.

9.4 Warm Starting

Warm starting a problem refers to providing a starting point x and additional information used by the solver to start the optimization process, for example, by providing information on which variables are suberbasic or active and thus hinting on the possible final active-set or providing a good initial guess for the final values of the Lagrange multipliers.
In order to warm start e04srf, it is necessary to
  1. (i)provide on the call to e04srf the initial guess x0;
  2. (ii)provide on the call to e04srf the initial guess for the Lagrange multipliers u. If nnzu>0 then the solver will access array u and so it must be provided. If nothing is known about the multipliers then u should be set to zero in the call to e04srf;
  3. (iii)store in the handle (under the label ‘SLACK’) the slack variables array. See Section 9.2;
  4. (iv)store in the handle (under the label ‘BASIS’ or ‘WARM START BASIS’) a valid basis vector of length nvar + stats(30). See Section 9.3;
  5. (v)request the solver to attempt a warm start by setting optional parameter SSQP Start Type=WARM.
If optional parameter SSQP Start Type=WARM but e04srf does not find the required information or it is inconsistent, then it will revert to a cold start.
Note: e04srf at exit (if the information is available) stores the slack variables and basis arrays into the handle under the labels ‘SLACK’ and ‘BASIS’. A next call to e04srf with the same handle along with SSQP Start Type=WARM, and the latest x and u, should trigger a warm start successfully. It will also notify the source of the warm starting information with a message similar to:
 Warm start information loaded successfully from handle.
       Handle [WARM START BASIS] data origin: solver
       Handle [SLACK] data origin           : user
Which indicates that the warm start information was successfully loaded. It also informs that you provided the slack variables array (using e04rxf with cmdstr=‘SLACK’) and that the basis information was provided by the solver itself, say, from a previous call to e04srf.

9.5 Internal Changes

Internal changes have been made to this routine as follows:
For details of all known issues which have been reported for the NAG Library please refer to the Known Issues.

10 Example

This example involves the minimization of the nonlinear function
f(x) = (x1+x2+x3)2 + 3x3 + 5x4 + cos (0.01x1) - 1 ,  
subject to the bounds
0 x3, 0 x4,  
to the nonlinear constraints
x12 + x22 + x3 = 2, x24 + x4 = 4,  
and the linear constraint
2x1 + 4x2 0 .  
The initial point is
x0 = (1,2,3,4)T ,  
and f(x0)=65. The optimal solution (to five significant figures) is
x* = (−0.070639,1.4124,0.0,0.019934) T , with ​ f(x*) = 1.9001 .  
Note: the variable x4 appears only linearly in the example. The linear property of the variable can be hinted to the solver using the routine e04rcf, see Section 3.1 and Section 3.1 in e04rcf for further details.

10.1 Program Text

Program Text (e04srfe.f90)

10.2 Program Results

Program Results (e04srfe.r)

11 Algorithmic Details

Here we summarise the main features of the SQP algorithm used in e04srf and introduce some terminology used in the description of the subroutine and its arguments. The SQP algorithm is fully described in Gill et al. (2002).

11.1 Constraints and Slack Variables

Let Q(x) be the mQ-dimensional vector made up of the mQ quadratic constraints. Furthermore, let m be the total number of components of the nonlinear constraints, g(x), the quadratic constraints, Q(x), and the linear constraints, Bx combined, that is, m=mg +mQ +mB . The upper and lower bounds on those terms define the general constraints of the problem. e04srf converts the general constraints to equalities by introducing a set of slack variables s = (s1,s2,,sm) T . For example, the linear constraint 5 2x1 + 3x2 is replaced by 2x1 + 3x2 - s1 = 0 together with the bounded slack 5 s1 . The minimization problem (1) can, therefore, be written in the equivalent form
minimize x,s f(x)  subject to ​ ( g(x) Q(x) Bx ) - s = 0 ,  l ( x s ) u . (2)
The general constraints become the equalities (g(x),Q(x))T - sN = 0 and Bx - sL = 0 , where sL and sN are the linear and nonlinear slacks.
Note: see Section 9.2 on how to access the slack variables stored in the handle.

11.2 Major Iterations

The basic structure of the SQP algorithm involves major and minor iterations. The major iterations generate a sequence of iterates {xk} that satisfy the linear constraints and converge to a point that satisfies the nonlinear constraints and the first-order conditions for optimality. At each iterate xk a QP subproblem is used to generate a search direction towards the next iterate xk+1 .
Without loss of generality and for ease of exposition, suppose that the model to optimize consists of linear and nonlinear constraints only, that is, mQ=0, then the constraints of the subproblem are formed from the linear constraints Bx - sL = 0 and the linearized constraint
g(xk) + g (xk) (x-xk) - sN = 0 , (3)
where g (xk) denotes the Jacobian matrix, whose elements are the first derivatives of g(x) evaluated at xk . The QP constraints, therefore, comprise the m linear constraints
g (xk) x - sN = -g(xk) + g (xk) xk , Bx -sN - sL = 0 , (4)
where x and s are bounded above and below by u and l as before. If the m × n matrix A and m -vector b are defined as
A = ( g (xk) B )  and  b = ( -g(xk) + g (xk) xk 0 ) , (5)
then the QP subproblem can be written as
minimize x,s Q (x,xk) = gkT (x-xk) + 12 (x-xk)T Hk (x-xk)  subject to ​ Ax - s = b ,  l ( x s ) u , (6)
where Q(x,xk) is a quadratic approximation to a modified Lagrangian function (see Gill et al. (2002)), g is the gradient of the Lagrangian, and the matrix Hk 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 nNLN rows and columns of the Hessian need to be approximated, where nNLN is the number of nonlinear variables, see Section 3.1 and Section 3.1 in e04rcf.

11.3 Minor Iterations

Solving the QP subproblem is itself an iterative procedure. Here, the iterations of the QP solver e04nqf form the minor iterations of the SQP method. e04nqf uses a reduced-Hessian active-set method implemented as a reduced-gradient method. At each minor iteration, the constraints Ax-s=b are partitioned into the form
BxB + SxS + NxN = b , (7)
where the basis matrix B is square and nonsingular, and the matrices S and N are the remaining columns of ( A -I ) . The vectors xB , xS and xN are the associated basic, superbasic and nonbasic variables respectively; they are a permutation of the elements of x and s . At a QP subproblem, the basic and superbasic variables will lie somewhere between their bounds, while the nonbasic variables will normally be equal to one of their bounds. At each iteration, xS is regarded as a set of independent variables that are free to move in any desired direction, namely one that will improve the value of the QP objective (or the sum of infeasibilities). The basic variables are then adjusted in order to ensure that (x,s) continues to satisfy Ax-s=b . The number of superbasic variables ( nS , say), therefore, indicates the number of degrees of freedom remaining after the constraints have been satisfied. In broad terms, nS is a measure of how nonlinear the problem is. In particular, nS will always be zero for LP problems.
If it appears that no improvement can be made with the current definition of B , S and N , a nonbasic variable is selected to be added to S , and the process is repeated with the value of nS increased by one. At all stages, if a basic or superbasic variable encounters one of its bounds, the variable is made nonbasic and the value of nS is decreased by one.
Associated with each of the m equality constraints Ax-s=b are the dual variables π. Similarly, each variable in (x,s) has an associated reduced gradient dj . The reduced gradients for the variables x are the quantities Q -ATπ , where Q is the gradient of the QP objective (6), and the reduced gradients for the slacks are the dual variables π . The QP subproblem is optimal if dj0 for all nonbasic variables at their lower bounds, dj0 for all nonbasic variables at their upper bounds, and dj=0 for other variables, including superbasics. In practice, an approximate QP solution (x^k,s^k,π^k) is found by relaxing these conditions.

11.4 The Merit Function

After a QP subproblem has been solved, new estimates of the solution are computed using a linesearch on the augmented Lagrangian merit function
M (x,s,π) = f(x) - πT (g(x)-sN) + 12 (g(x)-sN) T D (f(x)-sN) , (8)
where D is a diagonal matrix of penalty parameters (Dii0) , and π now refers to dual variables for the nonlinear constraints in (1). If (xk,sk,πk) denotes the current solution estimate and (x^k,s^k,π^k) denotes the QP solution, the linesearch determines a step αk (0<αk1) such that the new point
( xk+1 sk+1 πk+1 ) = ( xk sk πk ) + αk ( x^k - xk s^k - sk π^k - πk ) (9)
gives a sufficient decrease in the merit function M. When necessary, the penalties in D are increased by the minimum-norm perturbation that ensures descent for M (see Gill et al. (1992)). The value of sN is adjusted to minimize the merit function as a function of s before the solution of the QP subproblem (see Gill et al. (1986) and Eldersveld (1991)).

11.5 Treatment of Constraint Infeasibilities

e04srf makes explicit allowance for infeasible constraints. First, infeasible linear constraints are detected by solving the linear program
minimize x,v,w eT (v+w)  subject to ​ l ( x Bx - v + w ) u ,  v0 ,  w0 , (10)
where e is a vector of ones, and the nonlinear constraint bounds are temporarily excluded from l and u . This is equivalent to minimizing the sum of the general linear constraint violations subject to the bounds on x . (The sum is the 1 -norm of the linear constraint violations. In the linear programming literature, the approach is called elastic programming.)
The linear constraints are infeasible if the optimal solution of (10) has v0 or w0 . e04srf then terminates without computing the nonlinear functions.
Otherwise, all subsequent iterates satisfy the linear constraints. (Such a strategy allows linear constraints to be used to define a region in which the nonlinear functions can be safely evaluated.) e04srf proceeds to solve the nonlinear problem as given, using search directions obtained from the sequence of QP subproblems, see (6).
If a QP subproblem proves to be infeasible or unbounded (or if the dual variables π for the nonlinear constraints become large), e04srf enters ‘elastic’ mode and thereafter solves the problem
minimize x,v,w f(x) + γ eT (v+w)  subject to ​ l ( x g(x) - v + w Bx ) u ,  v0 ,  w0 , (11)
where γ is a non-negative argument (the elastic weight), and f(x) + γ eT (v+w) is called a composite objective (the 1 penalty function for the nonlinear constraints).
The value of γ may increase automatically by multiples of 10 if the optimal v and w continue to be nonzero. If γ is sufficiently large, this is equivalent to minimizing the sum of the nonlinear constraint violations subject to the linear constraints and bounds.

12 Optional Parameters

Several optional parameters in e04srf define choices in the problem specification or the algorithm logic. In order to reduce the number of formal arguments of e04srf these optional parameters have associated default values that are appropriate for most problems. Therefore, you need only specify those optional parameters whose values are to be different from their default values.
The remainder of this section can be skipped if you wish to use the default values for all optional parameters.
The optional parameters can be changed by calling e04zmf anytime between the initialization of the handle by e04raf and the call to the solver. Modification of the arguments during intermediate monitoring steps is not allowed. Once the solver finishes, the optional parameters can be altered again for the next solve.
The option values may be retrieved by e04znf.
The following is a list of the optional parameters available. A full description of each optional parameter is provided in Section 12.1.

12.1 Description of the Optional Parameters

For each option, we give a summary line, a description of the optional parameter and details of constraints.
The summary line contains:
All options accept the value DEFAULT to return single options to their default states.
Keywords and character values are case and white space insensitive.
Defaults
This special keyword may be used to reset all optional parameters to their default values. Any value given with this keyword will be ignored.
Infinite Bound SizerDefault =1020
This defines the ‘infinite’ bound bigbnd in the definition of the problem constraints. Any upper bound greater than or equal to bigbnd will be regarded as + (and similarly any lower bound less than or equal to -bigbnd will be regarded as -). Note that a modification of this optional parameter does not influence constraints which have already been defined; only the constraints formulated after the change will be affected.
Constraint: Infinite Bound Size1000.
Monitoring FileiDefault =−1
If i0, the unit number for the secondary (monitoring) output. If set to −1, no secondary output is provided. The following information is output to the unit:
Constraint: Monitoring File−1.
Monitoring LeveliDefault =4
This parameter sets the amount of information detail that will be printed by the solver to the secondary output. The meaning of the levels is the same as with Print Level.
Constraint: 0Monitoring Level5.
Print FileiDefault =advisory message unit number
If i0, the unit number for the primary output of the solver. If Print File=−1, the primary output is completely turned off independently of other settings. The default value is the advisory message unit number as defined by x04abf at the time of the optional parameters initialization, e.g., at the initialization of the handle. The following information is output to the unit:
Constraint: Print File−1.
Print LeveliDefault =2
This parameter defines how detailed information should be printed by the solver to the primary output.
i Output
0 No output from the solver
1 Only the final status and the objective value
2 Problem statistics, one line per iteration showing the progress of the solution with respect to the convergence measures, final status and statistics
3 As level 2 and additionally inner QP iteration log is reported every SSQP Print Frequency
4,5 As level 3 but further details of each minor and major iteration are presented
Constraint: 0Print Level5.
Print OptionsaDefault =YES
If Print Options=YES, a listing of optional parameters will be printed to the primary output.
Constraint: Print Options=YES or NO.
Print SolutionaDefault =NO
If Print Solution=YES, the final values of the solution vector are printed on the primary and secondary outputs.
Constraint: Print Solution=YES or NO.
SSQP Crash OptionaDefault =TRIPLE-PASS
Except on restarts SSQP Start Type=COLD, a crash procedure is used to select an initial basis B (see (7) in Section 11.3) from certain rows and columns of the constraint matrix ( A -I ) . The SSQP Crash Option determines which rows and columns of A (see (5) in Section 11.2) are eligible initially, and how many times crash is called. Columns of -I are used to pad the basis where necessary.
a Meaning
ONLY-SLACKS The initial basis contains only slack variables: B=I (see (7) in Section 11.3).
SINGLE-PASS crash is called once, looking for a triangular basis in all rows and columns of A.
DOUBLE-PASS 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).
TRIPLE-PASS 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.
If SSQP Crash OptionONLY-SLACKS, certain slacks on inequality rows are selected for the basis first, furthermore if SSQP Crash Option=DOUBLE-PASS or TRIPLE-PASS, 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 remaining unassigned rows, slack variables are inserted to complete the basis.
Note: each crash pass is counted as a minor iteration..
Constraint: SSQP Crash Option=ONLY-SLACKS, SINGLE-PASS, DOUBLE-PASS or TRIPLE-PASS.
SSQP Crash TolerancerDefault =0.1
The SSQP Crash Tolerance r allows the starting procedure crash to ignore certain small nonzeros in each column of A (see (5) in Section 11.2). If amax is the largest element in column j, other nonzeros aij in the column are ignored if |aij|amax×r. To be meaningful, r needs to be in the range 0r1.)
When r>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 (columns of -I). A feasible solution may be reached sooner on some problems.
When r=-1 (or r<0) the solver will automatically choose an appropriate tolerance.
Constraint: SSQP Crash Tolerance0 or SSQP Crash Tolerance=−1.0.
SSQP Estimate DerivativesaDefault =NO
This option specifies which nonlinear function gradients (objective and nonlinear constraints) are known analytically and will be supplied to the solver by you subroutines objgrd and congrd. This option has the following alternatives:
a Meaning
NO (default) Indicates to the solver that you are providing full derivative information, that is, all objective and constraint gradients are known.
ONLYOBJ Indicates that some or all components of the objective gradient will be missing but all constraint gradients are known.
ONLYCON Conversely to ONLYOBJ, this indicated that some or all of the constraint gradients will be missing but the objective gradient is known.
YES Indicates that the user is not providing derivatives.
The value SSQP Estimate Derivatives=NO should be used whenever possible. It is the most reliable and will usually be the most efficient.
If SSQP Estimate Derivatives=YES or ONLYOBJ, e04srf will estimate the missing components of the objective gradient, using finite differences. This may simplify the coding of subroutine objgrd. However, it could increase the total run-time substantially (since a special call to objfun is required for each missing element), and there is less assurance that an acceptable solution will be located. If the nonlinear variables are not well scaled, it may be necessary to specify a nonstandard difference interval changing SSQP Finite Diff Ctrl Interval or SSQP Finite Diff Interval.
Conversely, if SSQP Estimate Derivatives=YES or ONLYCON, e04srf will estimate missing elements of the Jacobian. For each column of the Jacobian, one call to confun is needed to estimate all missing elements in that column, if any.
At times, central differences are used rather than forward differences. Twice as many calls to objfun and confun are then needed. (This is not under your control.)
Note: if SSQP Estimate Derivatives=NO, ONLYCON or ONLYOBJ and the user does not provide one or more of the required gradient (or Jacobian) elements then these will be assumed to be constant, this might not be the desired behaviour. Under this circumstance, e04srf will inform the user with a message.
Constraint: SSQP Estimate Derivatives=NO, ONLYOBJ, ONLYCON or YES.
SSQP Finite Diff Ctrl IntervalrDefault 1.7×ε0.27
When SSQP Estimate Derivatives=YES, ONLYOBJ or ONLYCON, the central-difference interval r is used near an optimal solution to obtain more accurate (but more expensive) estimates of gradients. Twice as many function evaluations are required compared to forward differencing. The interval used for the j th variable is hj=r (1+|xj|) . The resulting derivative estimates should be accurate to O(r2) , unless the functions are badly scaled.
Constraint: 0<r<0.1.
SSQP Finite Diff IntervalrDefault 1.7×ε0.4
This alters the interval r used to estimate gradients by forward differences. It does so in the following circumstances:
In all cases, a derivative with respect to xj is estimated by perturbing that component of x to the value xj + r (1+|xj|) , and then evaluating the objective or constraint functions at the perturbed point. The resulting gradient estimates should be accurate to O(r) unless the functions are badly scaled. Judicious alteration of r may sometimes lead to greater accuracy.
Constraint: 0<r<0.1.
SSQP Function PrecisionrDefault 1.7×ε0.8
The relative function precision ε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 εr would be 1.0E−6 .
Ideally, the functions f(x) and gi(x) should have a magnitude of order 1. If all functions are substantially less than 1 in magnitude, ε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 εr would be 1.0E−10 .)
The default value of εr is appropriate for simple analytic functions.
In some cases the function values will be the result of extensive computation, possibly involving a costly iterative procedure that can provide few digits of precision. Specifying an appropriate SSQP 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.
Constraint: 0<εr<0.1.
SSQP HessianaDefault =AUTO
This option selects the method for storing and updating the approximate Hessian. (e04srf uses a quasi-Newton approximation to the Hessian of the Lagrangian. A BFGS update is applied after each major iteration.)
If SSQP 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 variables n is not too large (say, less than 75). In this case, the storage requirement is fixed and one can expect 1-step Q-superlinear convergence to the solution.
SSQP Hessian=LIMITED-MEMORY should be used on problems where n is very large. In this case, a limited-memory procedure is used to update a diagonal Hessian approximation Hr a limited number of times. Updates are accumulated as a list of vector pairs. They are discarded at regular intervals after Hr has been reset to their diagonal, see optional parameter SSQP Hessian Updates.
If SSQP Hessian=AUTO is specified, then the quasi-Newton method used is chosen based on a heuristic involving the number of nonlinear variables (see Section 3.1) in the problem.
Constraint: SSQP Hessian=AUTO, FULL-MEMORY or LIMITED-MEMORY.
SSQP Hessian UpdatesiDefault =10
If optional parameter SSQP Hessian=LIMITED-MEMORY is in effect 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 more vectors stored, the greater the cost of each QP iteration. 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 ).
Constraint: 0<SSQP Hessian Updates.
SSQP Iteration LimitiDefault = max(10000, 10 max(nvar,m) )
The value of i specifies the maximum number of minor iterations allowed (i.e., iterations of the simplex method or the QP algorithm), summed over all major iterations. (See also the description of the optional parameter SSQP Minor Iteration Limit.) For the default value, m is the number of linear, quadratic and nonlinear constraints defined in the model.
Constraint: 0<SSQP Iteration Limit.
SSQP Major Feasibility TolrDefault ε0.38
This tolerance, r , specifies how accurately the nonlinear constraints should be satisfied. The default value is appropriate when the linear and nonlinear constraints contain data to about that accuracy.
Let vmax be the maximum nonlinear constraint violation, normalized by the size of the solution, which is required to satisfy
vmax=maxi vi ​ ​/​ ​ x r , (12)
where vi is the violation of the i th nonlinear constraint, for i=1,2,,ncnln.
In the major iteration log (see Section 9.1), vmax appears as the quantity labelled ‘Feasible’. If some of the problem functions are known to be of low accuracy, a larger SSQP Major Feasibility Tol may be appropriate.
Constraint: 0<r0.1.
SSQP Major Iteration LimitiDefault = max(1000, 3 max(nvar,m) )
This is the maximum number of major iterations allowed. It is intended to guard against an excessive number of linearizations of the constraints. For the default value, m is the number of linear, quadratic and nonlinear constraints defined in the model.
Constraint: 0<SSQP Major Iteration Limit.
SSQP Major Optimality TolrDefault ε0.38
This tolerance, r , specifies the final accuracy of the dual variables. On successful termination, e04srf will have computed a solution (x,s,π) such that
cmax=maxj cj ​ ​/​ ​ π r , (13)
where cj is an estimate of the complementarity slackness for variable j , for j=1,2,,n+m. Here m represents the number of constraints. The values ci are computed from the final QP solution using the reduced gradients dj=gj - πT aj (where gj is the j th component of the objective gradient, aj is the associated column of the constraint matrix ( A -I ), and π is the set of QP dual variables):
cj = { - dj min( xj - lj ,1) if ​ dj0 ; -dj min( uj - xj ,1) if ​ dj<0 . ) (14)
In the major iteration log, cmax appears as the quantity labelled ‘Optimal’, See section Section 9.1.5.
Constraint: 0<r0.1.
SSQP Minor Feasibility TolrDefault ε0.38
e04srf tries to ensure that all variables eventually satisfy their upper and lower bounds to within this tolerance, r . This includes slack variables. Hence, general linear constraints should also be satisfied to within r .
Feasibility with respect to nonlinear constraints is judged by the optional parameter SSQP Major Feasibility Tol (not by SSQP Minor Feasibility Tol).
If the bounds and linear constraints cannot be satisfied to within r , the problem is declared infeasible. If rinfo(12) is quite small, it may be appropriate to raise r 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, except perhaps in presolve where the solver may require to evaluate confun, objgrd or congrd at a bound-feasible point. 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)=x1 + log(x2) , it is essential to place lower bounds on both variables. If r=1.0E−6 , the bounds x1 10−5 and x2 10−4 might be appropriate. (The log singularity is more serious. In general, keep x as far away from singularities as possible.)
If SSQP Scale Option=ALL or VARS AND LC, feasibility is defined in terms of the scaled problem (since it is then more likely to be meaningful).
In reality, e04srf uses r 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. e04srf then activates elastic mode for the rest of the solving process (with only the linearized nonlinear constraints defined to be elastic). For further details see Section 11.5.
Constraint: 0<r0.1.
SSQP Minor Iteration LimitiDefault =500
If the number of minor iterations for the optimality phase of the QP subproblem exceeds i , 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 i 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 (see Section 9.1) the flag t at the end of a line indicates that the corresponding QP was artificially terminated using the limit i .
Compare with the optional parameter SSQP Iteration Limit, which defines an independent absolute limit on the total number of minor iterations (summed over all QP subproblems).
Constraint: 0<SSQP Minor Iteration Limit.
SSQP Monitor FrequencyiDefault =0
This parameter specifies the frequency on which to call the monitor routine monit. If zero, the monitor routine will not be called.
Constraint: SSQP Monitor Frequency0.
SSQP Penalty ParameterrDefault =0.0
Defines the initial value of the penalty parameter and sets D=diag(r) in (8).
Constraint: 0.0SSQP Penalty Parameter.
SSQP Print FrequencyiDefault =1
Defines how frequently to print log information for the inner iterations. Only relevant if Print Level3 (with Print File−1) or Monitoring Level3 (with Monitoring File−1)
Constraint: 0SSQP Print Frequency.
SSQP Scale OptionaDefault =NONE
Three scale options are available as follows:
a Meaning
NONE No scaling. This is recommended if it is known that x and the constraint matrix never have very large elements (say, larger than 100).
VARS AND LC The 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 (1982)). This will sometimes improve the performance of the solution procedures.
ALL The constraints and variables are scaled by the iterative procedure. Also, a certain additional scaling is performed that may be helpful if the right-hand side b or the solution x is large. This takes into account columns of ( A -I ) that are fixed or have positive lower bounds or negative upper bounds.
Constraint: SSQP Scale Option=NONE, VARS AND LC or ALL.
SSQP Scale PrintaDefault =NO
Defines whether to print the row scales r(i) and column scales c(j). The scaled matrix coefficients are a¯ ij = aij c(j) / r(i) and the scaled bounds on the variables and slacks are l¯j = lj / c(j) , u¯j = uj / r(j-n) if j>n.
Constraint: SSQP Scale Print=YES or NO.
SSQP Scale TolerancerDefault =0.9
This parameter 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:
ρj=maxj |aij| / mini |aij| (aij0) .  
If maxj ρ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. The value of r should lie in the range 0<r<1.
Constraint: εSSQP Scale Tolerance.
SSQP Start TypeaDefault =COLD
Defines whether to perform a cold or warm start. If warm start data is not provided or is considered to have an unexpected size or content, then the solver will revert to perform a cold start on the problem. See Section 9.4 on how to correctly warm start a problem.
Constraint: SSQP Start Type=COLD or WARM.
SSQP Superbasics LimitiDefault =−1
This option 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 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 .
If i=−1 (the default), then e04srf automatically sets the limit to the number of nonlinear variables. For many problems, i may be considerably smaller than n. This will save storage if n is very large.
If the limit set is too low then the solver will terminate with ifail=27 error message.
Constraint: −1SSQP Superbasics Limit.
Stats TimeaDefault =NO
This parameter allows you to turn on the timings of various parts of the algorithm to give a better overview of where most of the time is spent. This might be helpful for a choice of different solving approaches. It is possible to choose between CPU and wall clock time. Choice YES is equivalent to WALL CLOCK.
Constraint: Stats Time=YES, NO, CPU or WALL CLOCK.
TaskaDefault =MINIMIZE
This parameter specifies the required direction of the optimization. If Task=FEASIBLE POINT, the objective function (if set) is ignored and the algorithm stops as soon as a feasible point is found. If no objective function is set, Task reverts to FEASIBLE POINT automatically.
Constraint: Task=MINIMIZE, MAXIMIZE or FEASIBLE POINT.
Time LimitrDefault =106
This parameter specifies a limit in seconds that the solver can use to solve one problem. If during the convergence check this limit is exceeded, the solver will terminate with ifail=23 error message.
Constraint: Time Limit>0.
Verify DerivativesaDefault =NO
This parameter specifies whether the routine should perform numerical checks on the consistency of the user-supplied gradient functions objgrd and congrd. If any discrepancies are found, ifail=26 is returned. It is recommended that such checks are enabled when first developing the formulation of the problem, however, the verification process results in a significant increase in the number of the function evaluations and thus it shouldn't be used in production code.
a Meaning
CHEAP Only a ‘cheap’ test will be performed.
YES Individual columns of the problem Jacobian, as well as individual gradient elements, will be checked (with a more reliable test). A key of the form OK or Bad? indicates whether or not each component appears to be correct.
NO Derivative checking is disabled.
Verify Derivatives=YES should be specified whenever a new objgrd or congrd is being developed.
Constraint: Verify Derivatives=CHEAP, YES or NO.