e04vh Method
関数リスト一覧 　　NagLibrary Namespaceへ　　ライブラリイントロダクション　　本ヘルプドキュメントのchm形式版

e04vh solves sparse nonlinear programming problems.

Syntax

C#
public static void e04vh(
	int start,
	int nf,
	int n,
	int nxname,
	int nfname,
	double objadd,
	int objrow,
	string prob,
	E04..::.E04VH_USRFUN usrfun,
	int[] iafun,
	int[] javar,
	int nea,
	double[] a,
	int[] igfun,
	int[] jgvar,
	int neg,
	double[] xlow,
	double[] xupp,
	string[] xnames,
	double[] flow,
	double[] fupp,
	string[] fnames,
	double[] x,
	int[] xstate,
	double[] xmul,
	double[] f,
	int[] fstate,
	double[] fmul,
	ref int ns,
	out int ninf,
	out double sinf,
	E04..::.e04vhOptions options,
	out int ifail
)
Visual Basic (Declaration)
Public Shared Sub e04vh ( _
	start As Integer, _
	nf As Integer, _
	n As Integer, _
	nxname As Integer, _
	nfname As Integer, _
	objadd As Double, _
	objrow As Integer, _
	prob As String, _
	usrfun As E04..::.E04VH_USRFUN, _
	iafun As Integer(), _
	javar As Integer(), _
	nea As Integer, _
	a As Double(), _
	igfun As Integer(), _
	jgvar As Integer(), _
	neg As Integer, _
	xlow As Double(), _
	xupp As Double(), _
	xnames As String(), _
	flow As Double(), _
	fupp As Double(), _
	fnames As String(), _
	x As Double(), _
	xstate As Integer(), _
	xmul As Double(), _
	f As Double(), _
	fstate As Integer(), _
	fmul As Double(), _
	ByRef ns As Integer, _
	<OutAttribute> ByRef ninf As Integer, _
	<OutAttribute> ByRef sinf As Double, _
	options As E04..::.e04vhOptions, _
	<OutAttribute> ByRef ifail As Integer _
)
Visual C++
public:
static void e04vh(
	int start, 
	int nf, 
	int n, 
	int nxname, 
	int nfname, 
	double objadd, 
	int objrow, 
	String^ prob, 
	E04..::.E04VH_USRFUN^ usrfun, 
	array<int>^ iafun, 
	array<int>^ javar, 
	int nea, 
	array<double>^ a, 
	array<int>^ igfun, 
	array<int>^ jgvar, 
	int neg, 
	array<double>^ xlow, 
	array<double>^ xupp, 
	array<String^>^ xnames, 
	array<double>^ flow, 
	array<double>^ fupp, 
	array<String^>^ fnames, 
	array<double>^ x, 
	array<int>^ xstate, 
	array<double>^ xmul, 
	array<double>^ f, 
	array<int>^ fstate, 
	array<double>^ fmul, 
	int% ns, 
	[OutAttribute] int% ninf, 
	[OutAttribute] double% sinf, 
	E04..::.e04vhOptions^ options, 
	[OutAttribute] int% ifail
)
F#
static member e04vh : 
        start:int * 
        nf:int * 
        n:int * 
        nxname:int * 
        nfname:int * 
        objadd:float * 
        objrow:int * 
        prob:string * 
        usrfun:E04..::.E04VH_USRFUN * 
        iafun:int[] * 
        javar:int[] * 
        nea:int * 
        a:float[] * 
        igfun:int[] * 
        jgvar:int[] * 
        neg:int * 
        xlow:float[] * 
        xupp:float[] * 
        xnames:string[] * 
        flow:float[] * 
        fupp:float[] * 
        fnames:string[] * 
        x:float[] * 
        xstate:int[] * 
        xmul:float[] * 
        f:float[] * 
        fstate:int[] * 
        fmul:float[] * 
        ns:int byref * 
        ninf:int byref * 
        sinf:float byref * 
        options:E04..::.e04vhOptions * 
        ifail:int byref -> unit

Parameters

start
Type: System..::.Int32
On entry: indicates how a starting point is to be obtained.
start=0
Requests that the Crash procedure be used, unless a Basis file is provided via optional parameters Old Basis File, Insert File or Load File.
start=1
Is the same as start=0 but is more meaningful when a Basis file is given.
start=2
Means that xstate and fstate define a valid starting point (probably from an earlier call, though not necessarily).
Constraint: start=0, 1 or 2.
nf
Type: System..::.Int32
On entry: nf, the number of problem functions in Fx , including the objective function (if any) and the linear and nonlinear constraints. Upper and lower bounds on x  can be defined using the parameters xlow and xupp and should not be included in F .
Constraint: nf>0.
n
Type: System..::.Int32
On entry: n , the number of variables.
Constraint: n>0.
nxname
Type: System..::.Int32
On entry: the number of names provided in the array xnames.
nxname=1
There are no names provided and generic names will be used in the output.
nxname=n
Names for all variables must be provided and will be used in the output.
Constraint: nxname=1​ or ​n.
nfname
Type: System..::.Int32
On entry: the number of names provided in the array fnames.
nfname=1
There are no names provided and generic names will be used in the output.
nfname=nf
Names for all functions must be provided and will be used in the output.
Constraint: nfname=1​ or ​nf.
Note: if nxname=1 then nfname must also be 1 (and vice versa). Similarly, if nxname=n then nfname must be nf (and vice versa).
objadd
Type: System..::.Double
On entry: is a constant that will be added to the objective row Fobj for printing purposes. Typically, objadd=0.0D+0.
objrow
Type: System..::.Int32
On entry: says which row of Fx is to act as the objective function. If there is no such row, set objrow=0 . Then e04vh will seek a feasible point such that lF ≤ Fx ≤ uF  and lx ≤ x ≤ ux .
Constraint: 1≤objrow≤nf​ or ​objrow=0 (or a feasible point problem).
prob
Type: System..::.String
On entry: is an 8-character name for the problem. prob is used in the printed solution and in some methods that output Basis files. A blank name may be used.
usrfun
Type: NagLibrary..::.E04..::.E04VH_USRFUN
usrfun must define the nonlinear portion fx  of the problem functions Fx=fx + Ax , along with its gradient elements Gij x= ∂ fix ∂ xj . (A dummy method is needed even if f≡0  and all functions are linear.)
In general, usrfun should return all function and gradient values on every entry except perhaps the last. This provides maximum reliability and corresponds to the default option setting, Derivative Option=1.
The elements of Gx  are stored in the array g[i-1], for i=1,2,…,leng, in the order specified by the input arrays igfun and jgvar.
In practice it is often convenient not to code gradients. e04vh is able to estimate them by finite differences, using a call to usrfun for each variable xj  for which some ∂ fix ∂ xj  needs to be estimated. However, this reduces the reliability of the optimization algorithm, and it can be very expensive if there are many such variables xj .
As a compromise, e04vh allows you to code as many gradients as you like. This option is implemented as follows. Just before usrfun is called, each element of the derivative array g is initialized to a specific value. On exit, any element retaining that value must be estimated by finite differences.
Some rules of thumb follow:
(i) for maximum reliability, compute all gradients;
(ii) if the gradients are expensive to compute, specify optional parameter Nonderivative Linesearch and use the value of the input parameter needg to avoid computing them on certain entries. (There is no need to compute gradients if needg=0 on entry to usrfun.);
(iii) if not all gradients are known, you must specify Derivative Option=0. You should still compute as many gradients as you can. (It often happens that some of them are constant or zero.);
(iv) again, if the known gradients are expensive, don't compute them if needg=0 on entry to usrfun;
(v) use the input parameter status to test for special actions on the first or last entries;
(vi) while usrfun is being developed, use the optional parameter Verify Level to check the computation of gradients that are supposedly known;
(vii) usrfun is not called until the linear constraints and bounds on x  are satisfied. This helps confine x  to regions where the functions fix  are likely to be defined. However, be aware of the optional parameter Minor Feasibility Tolerance if the functions have singularities on the constraint boundaries;
(viii) set status=-1  if some of the functions are undefined. The linesearch will shorten the step and try again;
(ix) set status≤-2  if you want e04vh to stop.

A delegate of type E04VH_USRFUN.

iafun
Type: array< System..::.Int32 >[]()[]
An array of size [lena]
Note: lena must satisfy the constraint: lena≥1
On entry: define the coordinates i,j  and values Aij  of the nonzero elements of the linear part A  of the function Fx=fx + Ax .
In particular, nea triples iafun[k-1],javar[k-1],a[k-1] define the row and column indices i=iafun[k-1] and j=javar[k-1] of the element Aij=a[k-1].
The coordinates may define the elements of A  in any order.
javar
Type: array< System..::.Int32 >[]()[]
An array of size [lena]
Note: lena must satisfy the constraint: lena≥1
On entry: define the coordinates i,j  and values Aij  of the nonzero elements of the linear part A  of the function Fx=fx + Ax .
In particular, nea triples iafun[k-1],javar[k-1],a[k-1] define the row and column indices i=iafun[k-1] and j=javar[k-1] of the element Aij=a[k-1].
The coordinates may define the elements of A  in any order.
nea
Type: System..::.Int32
On entry: is the number of nonzero entries in A  such that Fx=fx + Ax .
Constraint: 0≤nea≤lena.
a
Type: array< System..::.Double >[]()[]
An array of size [lena]
Note: lena must satisfy the constraint: lena≥1
On entry: define the coordinates i,j  and values Aij  of the nonzero elements of the linear part A  of the function Fx=fx + Ax .
In particular, nea triples iafun[k-1],javar[k-1],a[k-1] define the row and column indices i=iafun[k-1] and j=javar[k-1] of the element Aij=a[k-1].
The coordinates may define the elements of A  in any order.
igfun
Type: array< System..::.Int32 >[]()[]
An array of size [leng]
Note: leng must satisfy the constraint: leng≥1
On entry: define the coordinates i,j  of the nonzero elements of G , the nonlinear part of the derivative Jx=Gx + A  of the function Fx=fx + Ax . e04vj may be used to define these two arrays.
The coordinates can define the elements of G  in any order. However, usrfun must define the actual elements of g in exactly the same order as defined by the coordinates igfun,jgvar .
jgvar
Type: array< System..::.Int32 >[]()[]
An array of size [leng]
Note: leng must satisfy the constraint: leng≥1
On entry: define the coordinates i,j  of the nonzero elements of G , the nonlinear part of the derivative Jx=Gx + A  of the function Fx=fx + Ax . e04vj may be used to define these two arrays.
The coordinates can define the elements of G  in any order. However, usrfun must define the actual elements of g in exactly the same order as defined by the coordinates igfun,jgvar .
neg
Type: System..::.Int32
On entry: the number of nonzero entries in G .
Constraint: 0≤neg≤leng.
xlow
Type: array< System..::.Double >[]()[]
An array of size [n]
On entry: contain the lower and upper bounds lx  and ux  on the variables x .
To specify a nonexistent lower bound lx j = -∞ , set xlow[j-1]≤-infbnd , where infbnd  is the optional parameter Infinite Bound Size. To specify a nonexistent upper bound ux j = ∞ , set xupp[j-1]≥infbnd.
To fix the jth variable at xj=β , where β<infbnd, set xlow[j-1]=xupp[j-1]=β .
Constraint: xlow[i-1]≤xupp[i-1], for i=1,2,…,n.
xupp
Type: array< System..::.Double >[]()[]
An array of size [n]
On entry: contain the lower and upper bounds lx  and ux  on the variables x .
To specify a nonexistent lower bound lx j = -∞ , set xlow[j-1]≤-infbnd , where infbnd  is the optional parameter Infinite Bound Size. To specify a nonexistent upper bound ux j = ∞ , set xupp[j-1]≥infbnd.
To fix the jth variable at xj=β , where β<infbnd, set xlow[j-1]=xupp[j-1]=β .
Constraint: xlow[i-1]≤xupp[i-1], for i=1,2,…,n.
xnames
Type: array< System..::.String >[]()[]
An array of size [nxname]
On entry: the optional names for the variables.
If nxname=1 , xnames is not referenced and default names will be used for output.
If nxname=n , xnames[j-1]  should contain the 8-character name of the j th variable.
flow
Type: array< System..::.Double >[]()[]
An array of size [nf]
On entry: contain the lower and upper bounds lF  and uF  on Fx .
To specify a nonexistent lower bound lF i =-∞ , set flow[i-1]≤-infbnd. For a nonexistent upper bound uF i =∞ , set fupp[i-1]≥infbnd .
To make the i th constraint an equality at Fi=β , where β < infbnd , set flow[i-1]=fupp[i-1]=β.
Constraint: flow[i-1]≤fupp[i-1], for i=1,2,…,n.
fupp
Type: array< System..::.Double >[]()[]
An array of size [nf]
On entry: contain the lower and upper bounds lF  and uF  on Fx .
To specify a nonexistent lower bound lF i =-∞ , set flow[i-1]≤-infbnd. For a nonexistent upper bound uF i =∞ , set fupp[i-1]≥infbnd .
To make the i th constraint an equality at Fi=β , where β < infbnd , set flow[i-1]=fupp[i-1]=β.
Constraint: flow[i-1]≤fupp[i-1], for i=1,2,…,n.
fnames
Type: array< System..::.String >[]()[]
An array of size [nfname]
On entry: the optional names for the problem functions.
If nfname=1 , fnames is not referenced and default names will be used for output.
If nfname=nf , fnames[i-1]  should contain the 8-character name of the i th row of F .
x
Type: array< System..::.Double >[]()[]
An array of size [n]
On entry: an initial estimate of the variables x . See the following description of xstate.
On exit: the final values of the variable x .
xstate
Type: array< System..::.Int32 >[]()[]
An array of size [n]
On entry: the initial state for each variable x .
If start=0 or 1 and no basis information is provided (the optional parameters Old Basis File, Insert File and Load File are all set to 0; the default) x and xstate must be defined.
If nothing special is known about the problem, or if there is no wish to provide special information, you may set x[j-1]=0.0, xstate[j-1]=0, for all j=1,2,…,n . If you set x[j-1]=xlow[j-1]  set xstate[j-1]=4 ; if you set x[j-1]=xupp[j-1]  then set xstate[j-1]=5 . In this case a Crash procedure is used to select an initial basis.
If start=0 or 1 and basis information is provided (at least one of the optional parameters Old Basis File, Insert File and Load File is nonzero) x and xstate need not be set.
If start=2 (Warm Start), x and xstate must be set (probably from a previous call). In this case xstate[j-1] must be 0, 1, 2 ​ or ​ 3, for j=1,2,…,n.
On exit: the final state of the variables.
xstate[j-1] State of variable j Usual value of x[j-1]
0 nonbasic xlow[j-1]
1 nonbasic xupp[j-1]
2 superbasic Between xlow[j-1] and xupp[j-1]
3 basic Between xlow[j-1] and xupp[j-1]
Basic and superbasic variables may be outside their bounds by as much as the optional parameter Minor Feasibility Tolerance. Note that if scaling is specified, 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 value of Primal infeasibility output to the unit number associated with the optional parameter Print File.
Very occasionally some nonbasic variables may be outside their bounds by as much as the optional parameter Minor Feasibility Tolerance, and there may be some nonbasics for which x[j-1] lies strictly between its bounds.
If ninf>0 , some basic and superbasic variables may be outside their bounds by an arbitrary amount (bounded by sinf if scaling was not used).
Constraint: 0≤xstate[j-1]≤5, for ​j=1,2,…,n.
xmul
Type: array< System..::.Double >[]()[]
An array of size [n]
On exit: the vector of the dual variables (Lagrange multipliers) for the simple bounds lx ≤ x ≤ ux .
f
Type: array< System..::.Double >[]()[]
An array of size [nf]
On entry: an initial value for the problem functions F . See the following description of fstate.
On exit: the final values for the problem functions F  (the values F  at the final point x).
fstate
Type: array< System..::.Int32 >[]()[]
An array of size [nf]
On entry: the initial state for the problem functions F.
If start=0 or 1 and no basis information is provided (the optional parameters Old Basis File, Insert File and Load File are all set to 0; the default, f and fstate must be defined.
If nothing special is known about the problem, or if there is no wish to provide special information, you may set f[i-1]=0.0, fstate[i-1]=0, for all i=1,2,…,nf. Less trivially, to say that the optimal value of function f[i-1]  will probably be equal to one of its bounds, set f[i-1]=flow[i-1] and fstate[i-1]=4 or f[i-1]=fupp[i-1] and fstate[i-1]=5 as appropriate. In this case a Crash procedure is used to select an initial basis.
If start=0 or 1 and basis information is provided (at least one of the optional parameters Old Basis File, Insert File and Load File is nonzero), f and fstate need not be set.
If start=2 (Warm Start), f and fstate must be set (probably from a previous call). In this case fstate[i-1]  must be 0, 1, 2 or 3, for i=1,2,…,nf .
On exit: the final state of the variables. The elements of fstate have the following meaning:
fstate[i-1] State of the corresponding
slack variable		 Usual value of f[i-1]
0 nonbasic flow[i-1]
1 nonbasic fupp[i-1]
2 superbasic Between flow[i-1] and fupp[i-1]
3 basic Between flow[i-1] and fupp[i-1]
Basic and superbasic slack variables may lead to the corresponding functions being outside their bounds by as much as the optional parameter Minor Feasibility Tolerance.
Very occasionally some functions may be outside their bounds by as much as the optional parameter Minor Feasibility Tolerance, and there may be some nonbasics for which f[i-1]  lies strictly between its bounds.
If ninf>0 , some basic and superbasic variables may be outside their bounds by an arbitrary amount (bounded by sinf if scaling was not used).
Constraint: 0≤fstate[i-1]≤5, for ​i=1,2,…,nf.
fmul
Type: array< System..::.Double >[]()[]
An array of size [nf]
On entry: an estimate of γ , the vector of Lagrange multipliers (shadow prices) for the constraints lF ≤ Fx ≤ uF . All nf components must be defined. If nothing is known about γ , set fmul[i-1]=0.0 , for i=1,2,…,nf. For warm start use the values from a previous call.
On exit: the vector of the dual variables (Lagrange multipliers) for the general constraints lF ≤ Fx ≤ uF  
ns
Type: System..::.Int32 %
On entry: the number of superbasic variables. ns need not be specified for cold starts, but should retain its value from a previous call when warm start is used.
On exit: the final number of superbasic variables.
ninf
Type: System..::.Int32 %
On exit: the number of the infeasibilities of constraints that lie outside one of their bounds by more than the optional parameter Minor Feasibility Tolerance before the solution is unscaled.
sinf
Type: System..::.Double %
On exit: the sum of the infeasibilities of constraints that lie outside one of their bounds by more than the optional parameter Minor Feasibility Tolerance before the solution is unscaled.
If any linear constraints are infeasible, x minimizes the sum of the infeasibilities of the linear constraints subject to the upper and lower bounds being satisfied. In this case ninf gives the number of variables and linear constraints lying outside their upper or lower bounds. The nonlinear constraints are not evaluated.
Otherwise, x minimizes the sum of infeasibilities of the nonlinear constraints subject to the linear constraints and upper and lower bounds being satisfied. In this case ninf gives the number of components of Fx  lying outside their bounds by more than the optional parameter Minor Feasibility Tolerance. Again this is before the solution is unscaled.
options
Type: NagLibrary..::.E04..::.e04vhOptions
An object of type E04.e04vhOptions. Used to configure optional parameters to this method.
ifail
Type: System..::.Int32 %
On exit: ifail=0 unless the method detects an error (see [Error Indicators and Warnings]).
e04vh returns with ifail=0 if the iterates have converged to a point x that satisfies the first-order Kuhn–Tucker (see [Minor Iteration Log]) conditions to the accuracy requested by the optional parameter Major Optimality Tolerance, i.e., the projected gradient and active constraint residuals are negligible at x.
You should check whether the following four conditions are satisfied:
(i) the final value of rgNorm (see [Minor Iteration Log]) is significantly less than that at the starting point;
(ii) during the final major iterations, the values of Step and Minors (see [Major Iteration Log]) are both one;
(iii) the last few values of both rgNorm and SumInf (see [Minor Iteration Log]) become small at a fast linear rate; and
(iv) condHz (see [Major Iteration Log]) is small.
If all these conditions hold, x is almost certainly a local minimum of (1).
One caution 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 in the Print file (see [Description of Monitoring Information]) refers to the largest bound infeasibility and which variable is involved. If it is too large, consider restarting with a smaller Minor Feasibility Tolerance (say 10 times smaller) and perhaps Scale Option=0.
Similarly, Max Dual infeas in the Print file indicates which variable is most likely to be at a nonoptimal value. Broadly speaking, if
Max Dual infeas / Max pi = 10-d ,
then the objective function would probably change in the d th significant digit if optimization could be continued. If d  seems too large, consider restarting with a smaller Major Optimality Tolerance.
Finally, Nonlinear constraint violn in the Print file shows the maximum infeasibility for nonlinear rows. If it seems too large, consider restarting with a smaller Major Feasibility Tolerance.

Description

e04vh is designed to minimize a linear or nonlinear function subject to bounds on the variables and sparse linear or nonlinear constraints. It is suitable for large-scale linear and quadratic programming and for linearly constrained optimization, as well as for general nonlinear programs of the form
minimizex f0x   subject to ​ l ≤ x fx ALx ≤ u , (1)
where x  is an n -vector of variables, l  and u  are constant lower and upper bounds, f0 x  is a smooth scalar objective function, AL  is a sparse matrix, and fx  is a vector of smooth nonlinear constraint functions fix . The optional parameter Maximize specifies that f0x  should be maximized instead of minimized.
Ideally, the first derivatives (gradients) of f0x  and fix  should be known and coded by you. If only some of the gradients are known, e04vh estimates the missing ones by finite differences.
If f0x  is linear and fx  is absent, (1) is a linear program (LP) and e04vh applies the primal simplex method (see Dantzig (1963)). Sparse basis factors are maintained by LUSOL (see Gill et al. (1987)) as in MINOS (see Murtagh and Saunders (1995)).
If only the objective is nonlinear, the problem is linearly constrained (LC) and tends to solve more easily than the general case with nonlinear constraints (NC). For both nonlinear cases, e04vh applies a sparse sequential quadratic programming (SQP) method (see Gill et al. (2002)), using limited-memory quasi-Newton approximations to the Hessian of the Lagrangian. The merit function for step-length control is an augmented Lagrangian, as in the dense SQP solver e04wd (see Gill et al. (1986) and Gill et al. (1992)).
e04vh is suitable for nonlinear problems with thousands of constraints and variables, and is most efficient if only some of the variables enter nonlinearly, or there are relatively few degrees of freedom at a solution (i.e., many constraints are active). However, there is no limit on the number of degrees of freedom.
e04vh allows linear and nonlinear constraints and variables to be entered in an arbitrary order, and uses one method to define all the nonlinear functions.
The optimization problem is assumed to be in the form
minimizex Fobj x   subject to ​ lx ≤ x ≤ ux ,   lF ≤ Fx ≤ uF , (2)
where the upper and lower bounds are constant, Fx  is a vector of smooth linear and nonlinear constraint functions Fix , and Fobj x  is one of the components of F  to be minimized, as specified by the input parameter objrow. e04vh reorders the variables and constraints so that the problem is in the form (1).
Upper and lower bounds are specified for all variables and functions. The j th constraint may be defined as an equality by setting lj = uj . If certain bounds are not present, the associated elements of l or u should be set to special values that are treated as -∞ or +∞. Free variables and free constraints (‘free rows’) have both bounds infinite.
In general, the components of F  are structured in the sense that they are formed from sums of linear and nonlinear functions of just some of the variables. This structure can be exploited by e04vh.
In many cases, the vector Fx  is a sum of linear and nonlinear functions. e04vh allows these terms to be specified separately, so that the linear part is defined just once by the input arguments iafun, javar and a. Only the nonlinear part is recomputed at each x .
Suppose that each component of Fx  is of the form
Fix = fix + ∑ j=1 n Aij xj ,
where fix  is a nonlinear function (possibly zero) and the elements Aij  are constant. The nf  by n  Jacobian of Fx  is the sum of two sparse matrices of the same size: F′x = Gx + A , where Gx = f′x  and A  is the matrix with elements Aij . The two matrices must be non-overlapping in the sense that each element of the Jacobian F′x = Gx + A  comes from Gx  or A , but not both. The element cannot be split between Gx  and A .
For example, the function
Fx = 3x1 + ex2 x4 + x22 + 4x4 - x3 + x5 x2 + x32 + sin⁡x4 - 3x5 x1 - x3
can be written as
Fx = fx + Ax = ex2 x4 + x22 + 4x4 x32 + sin⁡x4 0 + 3x1 - x3 + x5 x2 - 3x5 x1 - x3 ,
in which case
F′x = 3 ex2 x4 + 2x2 -1 ex2 + 4 -1 0 1 2x3 cos⁡x4 -3 1 0 -1 0 -0
can be written as F′x = f′x + A = Gx + A , where
Gx = 0 ex2 x4 + 2x2 0 ex2 + 4 0 0 0 2x3 cos⁡x4 0 0 0 0 0 0 ,   A = 3 0 -1 0 -1 0 1 -0 0 -3 1 0 -1 0 -0 .
Note:  the element ex2+4  of F′x  appears in Gx  and is not split between Gx  and A  although it contains a linear term.
The nonzero elements of A  and G  are provided to e04vh in coordinate form. The elements of A  are entered as triples i,j,Aij  in the arrays iafun, javar and a. The sparsity pattern G  is entered as pairs i,j  in the arrays igfun and jgvar. The corresponding entries Gij  (any that are known) are assigned to appropriate array elements gk  in usrfun.
The elements of A  and G  may be stored in any order. Duplicate entries are ignored. igfun and jgvar may be defined automatically by method e04vj when Derivative Option=0 is specified and usrfun does not provide any gradients.
Throughout this document the symbol ε is used to represent the machine precision (see x02aj).
e04vh is based on SNOPTA, which is part of the SNOPT package described in Gill et al. (2005b).

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 (2005a) Users' guide for SQOPT 7: a Fortran package for large-scale linear and quadratic programming Report NA 05-1 Department of Mathematics, University of California, San Diego ftp://www.cam.ucsd.edu/pub/peg/reports/sqdoc7.pdf
Gill P E, Murray W and Saunders M A (2005b) Users' guide for SNOPT 7.1: a Fortran package for large-scale linear nonlinear programming Report NA 05-2 Department of Mathematics, University of California, San Diego ftp://www.cam.ucsd.edu/pub/peg/reports/sndoc7.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. Program.
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

Error Indicators and Warnings

Note: e04vh may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the method:
Some error messages may refer to parameters that are dropped from this interface (lena, leng, cw, lencw, iw, leniw, rw, lenrw, cuser, iuser, ruser) In these cases, an error in another parameter has usually caused an incorrect value to be inferred.
ifail=1
The initialization method (e04vg not in this release) has not been called or at least one of lencw, leniw and lenrw is less than 600.
ifail=2
An input parameter is invalid. The output message provides more details of the invalid argument.
ifail=3
Requested accuracy could not be achieved.
A feasible solution has been found, but the requested accuracy in the dual infeasibilities could not be achieved. An abnormal termination has occurred, but e04vh is within 10-2  of satisfying the Major Optimality Tolerance. Check that the Major Optimality Tolerance is not too small.
ifail=4
The problem appears to be infeasible.
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 (6) and (7) in [Major Iterations]), there is apparently no point that satisfies the bounds on x  and s . Violations as small as the Minor Feasibility Tolerance are ignored, but at least one component of x  or s  violates a bound by more than the tolerance.
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, e04vh 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), e04vh enters so-called ‘nonlinear elastic’ mode. The subproblem includes the original QP objective and the sum of the infeasibilities – suitably weighted using the optional parameter Elastic Weight. 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 eventually will be obtained for the perturbed constraints, and optimization can continue on the subproblem. If the nonlinear problem has no feasible solution, e04vh will tend to determine a ‘good’ infeasible point if the elastic weight is sufficiently large. (If the elastic weight were infinite, e04vh would locally minimize the nonlinear constraint violations subject to the linear constraints and bounds.)
Unfortunately, even though e04vh 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=5
The problem appears to be unbounded (or badly scaled).
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 Scale Option.
For nonlinear problems, e04vh 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 (as judged by the unbounded parameters (see [Major Iteration Log])), the problem is terminated and declared unbounded. To avoid large function values, it may be necessary to impose bounds on some of the variables in order to keep them away from singularities in the nonlinear functions.
The message may indicate an abnormal termination while enforcing the limit on the constraint violations. This exit implies that the objective is not bounded below in the feasible region defined by expanding the bounds by the value of the Violation Limit.
ifail=6
Iteration limit reached.
Either the Iterations Limit or the Major Iterations Limit was exceeded before the required solution could be found. Check the iteration log to be sure that progress was being made. If so, and if you caused a basis file to be saved by using the optional parameter New Basis File, consider restarting the run using the optional parameter Old Basis File to see whether further progress can be made. If you have no basis file available, you might rerun the problem after increasing the optional parameters Minor Iterations Limit and/or Major Iterations Limit.
If none of the above limits have been reached, 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 can't continue as the number of superbasic variables has already reached the limit specified by the optional parameter Superbasics Limit. In general, raise the Superbasics Limit  s  by a reasonable amount, bearing in mind the storage needed for the reduced Hessian.
ifail=7
Numerical difficulties have been encountered and no further progress can be made.
Several circumstances could lead to this exit.
  1. usrfun could be returning accurate function values but inaccurate gradients (or vice versa). This is the most likely cause. Study the comments given for ifail=8, and do your best to ensure that the coding is correct.
  2. The function and gradient values could be consistent, but their precision could be too low. For example, accidental use of a real data type when real was intended would lead to a relative function precision of about 10-6  instead of something like 10-15 . The default Major Optimality Tolerance 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. 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
    • Function Precision  t
    • Major Optimality Tolerance  t
    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. 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 Scale Option=1 if scaling has not yet been used.
    For certain highly structured basis matrices (notably those with band structure), a systematic growth may occur in the factor U . Consult the description of Umax and Growth in [Basis Factorization Statistics] and set the LU Factor Tolerance to 2.0 (or possibly even smaller, but not less than 1.0).
  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 must mean that the problem is badly scaled, or the LU Factor Tolerance is too much larger than 1.0. This is highly unlikely to occur.
ifail=8
Derivative appears to be incorrect.
A check has been made on some elements of the Jacobian as returned in the parameter g of usrfun. At least one value 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 e04vh 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=7.
Check the function and Jacobian computation very carefully in usrfun. 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 Level=0 to prevent individual elements from being checked. However, the optimization procedure may have difficulty.
ifail=9
Undefined user-supplied function.
You have indicated that the problem functions are undefined by assigning the value status=-1  on exit from usrfun. e04vh attempts to evaluate the problem functions closer to a point at which the functions are already known to be defined. This exit occurs if e04vh is unable to find a point at which the functions are defined. This will occur in the case of:
– undefined functions with no recovery possible;
– undefined functions at the first point;
– undefined functions at the first feasible point; or
– undefined functions when checking derivatives.
ifail=10
User requested termination.
You have indicated the wish to terminate solution of the current problem by setting status to a value < -1  on exit from usrfun.
ifail=11
Internal memory allocation failed when attempting to obtain the required workspace. Please contact NAG.
ifail=12
Internal memory allocation was insufficient. Please contact NAG.
ifail=13
An error has occurred in the basis package, perhaps indicating incorrect setup of arrays. Set the optional parameter Print File and examine the output carefully for further information.
ifail=14
An unexpected error has occurred. Set the optional parameter Print File and examine the output carefully for further information.
ifail=-1000
The array lengths are not the same for arrays
ifail=-8000
Negative dimension for array value
ifail=-6000
Invalid Parameters value

Accuracy

Further Comments

This section describes the final output produced by e04vh. Intermediate and other output are given in [Description of Monitoring Information].

The Final Output

If Print File > 0 , the final output, including a listing of the status of every variable and constraint will be sent to the channel numbers associated with Print File. The following describes the output for each constraint (row) and variable (column).

The ROWS section

General linear constraints take the form l≤ALx≤u . The i th constraint is therefore of the form
α≤ νix≤ β ,
where νi  is the ith row of AL .
Internally, the constraints take the form ALx-s=0 , where s is the set of slack variables (which satisfy the bounds l≤s≤u ). For the ith row it is the slack variable si  that is directly available and it is sometimes convenient to refer to its state. Nonlinear constraints α≤fix+νix≤β  are treated similarly, except that the row activity and degree of infeasibility are computed directly from fix + νix , rather than si .
A full stop (.) is printed for any numerical value that is exactly zero.
Label Description
Number is the value of n+i. (This is used internally to refer to si in the intermediate output.)
Row gives the name of the ith row.
State the state of the ith row relative to the bounds α  and β . The various states possible are as follows:
LL the row is at its lower limit, α .
UL the row is at its upper limit, β .
EQ the limits are the same ( α=β ).
FR si  is nonbasic and currently zero, even though it is free to take any value between its bounds α  and β .
BS si  is basic.
SBS si  is superbasic.
A key is sometimes printed before State. Note that unless the optional parameter Scale Option=0 is specified, the tests for assigning a key are applied to the variables of the scaled problem.
A Alternative optimum possible. The variable is nonbasic, but its reduced gradient is essentially zero. This means that if the variable were allowed to start moving away from its bound, there would be no change in the value of the objective function. The values of the other free variables might change, giving a genuine alternative solution. However, if there are any degenerate variables (labelled D), the actual change might prove to be zero, since one of them could encounter a bound immediately. In either case, the values of the Lagrange multipliers might also change.
D Degenerate. The variable is basic or superbasic, but it is equal (or very close) to one of its bounds.
I Infeasible. The variable is basic or superbasic and is currently violating one of its bounds by more than the value of the Feasibility Tolerance.
N Not precisely optimal. The variable is nonbasic or superbasic. If the value of the reduced gradient for the variable exceeds the value of the optional parameter Major Optimality Tolerance, the solution would not be declared optimal because the reduced gradient for the variable would not be considered negligible.
Activity is the value of νix  (or fix + νix  for nonlinear rows) at the final iterate.
Slack Activity is the value by which the row differs from its nearest bound. (For the free row (if any), it is set to Activity.)
Lower Limit is α, the lower bound on the row.
Upper Limit is β, the upper bound on the row.
Dual Activity is the value of the dual variable πi (the Lagrange multiplier for the i th constraint). The full vector π  always satisfies BTπ = gB , where B  is the current basis matrix and gB  contains the associated gradients for the current objective function. For FP problems, πi is set to zero.
i gives the index i of the ith row.

The COLUMNS section

Let the j th component of x  be the variable xj  and assume that it satisfies the bounds α ≤ xj ≤ β . A fullstop (.) is printed for any numerical value that is exactly zero.
Label Description
Number is the column number j. (This is used internally to refer to xj in the intermediate output.)
Column gives the name of xj.
State the state of xj  relative to the bounds α  and β . The various states possible are as follows:
LL xj  is nonbasic at its lower limit, α .
UL xj  is nonbasic at its upper limit, β .
EQ xj  is nonbasic and fixed at the value α = β .
FR xj  is nonbasic at some value strictly between its bounds: α< xj< β .
BS xj  is basic. Usually α< xj< β .
SBS xj  is superbasic. Usually α< xj< β .
A key is sometimes printed before State. Note that unless the optional parameter Scale Option=0 is specified, the tests for assigning a key are applied to the variables of the scaled problem.
A Alternative optimum possible. The variable is nonbasic, but its reduced gradient is essentially zero. This means that if the variable were allowed to start moving away from its bound, there would be no change in the value of the objective function. The values of the other free variables might change, giving a genuine alternative solution. However, if there are any degenerate variables (labelled D), the actual change might prove to be zero, since one of them could encounter a bound immediately. In either case, the values of the Lagrange multipliers might also change.
D Degenerate. The variable is basic or superbasic, but it is equal (or very close) to one of its bounds.
I Infeasible. The variable is basic or superbasic and is currently violating one of its bounds by more than the value of the Feasibility Tolerance.
N Not precisely optimal. The variable is nonbasic or superbasic. If the value of the reduced gradient for the variable exceeds the value of the optional parameter Major Optimality Tolerance, the solution would not be declared optimal because the reduced gradient for the variable would not be considered negligible.
Activity is the value of xj at the final iterate.
Obj Gradient is the value of gj at the final iterate. For FP problems, gj is set to zero.
Lower Limit is the lower bound specified for the variable. None indicates that xlow[j-1]≤-infbnd.
Upper Limit is the upper bound specified for the variable. None indicates that xupp[j-1]≥infbnd.
Reduced Gradnt is the value of the reduced gradient dj = gj - πT aj  where aj  is the j th column of the constraint matrix. For FP problems, dj is set to zero.
m + j is the value of m+j.
Note that movement off a constraint (as opposed to a variable moving away from its bound) can be interpreted as allowing the entry in the Slack Activity column to become positive.
Numerical values are output with a fixed number of digits; they are not guaranteed to be accurate to this precision.

Example

This example is a reformulation of Problem 74 from Hock and Schittkowski (1981) and involves the minimization of the nonlinear function
fx=10-6x33+23×10-6x43+3x3+2x4
subject to the bounds
-0.55≤x1≤ 0.55, -0.55≤x2≤ 0.55, - 0.00≤x3≤ 1200, - 0.00≤x4≤ 1200,
to the nonlinear constraints
1000 sin -x1 -0.25 + 1000 sin -x2 -0.25 - x3 = -894.8, 1000sinx1-0.25+1000sinx1-x2-0.25-x4 = -894.8, 1000sinx2-0.25+1000sinx2-x1-0.25 = -1294.8,
and to the linear constraints
-x1+x2≥-0.55, -x1-x2≥-0.55.
The initial point, which is infeasible, is
x0 = 0, 0, 0, 0 T ,
and fx0=0.
The optimal solution (to five figures) is
x*=0.11887,-0.39623,679.94,1026.0T,
and fx*=5126.4. All the nonlinear constraints are active at the solution.
The example in the document for e04vj solves the above problem. It first calls e04vj to determine the sparsity pattern before calling e04vh.
The example in the document for (e04vk not in this release) solves the above problem using some of the optional parameters described in [Optional Parameters].
The formulation of the problem combines the constraints and the objective into a single vector ( F ) which is split into linear part ( ALx ) and a nonlinear part ( f ). For example we could write
F = 1000 sin -x1 - 0.25 + 1000 sin -x2 - 0.25 - x3 1000 sin x1 - 0.25 + 1000 sin x1 - x2 - 0.25 - x4 1000 sin x2 - 0.25 + 1000 sin x2 - x1 - 0.25 -x1 + x2 -x1 - x2 10-6 x33 + 23 × 10-6 x43 + 3x3 + 2x4 = f + AL x
where
f = 1000 sin -x1 - 0.25 + 1000 sin -x2 - 0.25 1000 sin x1 - 0.25 + 1000 sin x1 - x2 - 0.25 1000 sin x2 - 0.25 + 1000 sin x2 - x1 - 0.25 0 0 10-6 x33 + 23 × 10-6 x43
and
AL = 0 0 -1 0 0 0 0 -1 0 0 0 0 -1 1 0 0 1 -1 0 0 0 0 3 2 .
The nonzero elements of AL  need to be stored in the triples iafun[k-1],javar[k-1],a[k-1]  in any order. For example
k 1 2 3 4 5 6 7 8
iafun[k-1] 1 2 4 4 5 5 6 6
javar[k-1] 3 4 1 2 1 2 3 4
a[k-1] -1 -1 -1 1 1 -1 3 2
The nonlinear functions f  and the Jacobian need to be supplied in usrfun. Please note that there is no need to assign any value to f4  or f5  as there is no nonlinear part in F4  or F5 .
The nonzero entries of the Jacobian of f  are
∂f1 ∂x1 = -1000 cos -x1 - 0.25 ∂f1 ∂x2 = -1000 cos -x2 - 0.25 ∂f2 ∂x1 = 1000 cos x1 - 0.25 + 1000 cos x1 - x2 - 0.25 ∂f2 ∂x2 = -1000 cos x1 - x2 - 0.25 ∂f3 ∂x1 = -1000 cos x2 - x1 - 0.25 ∂f3 ∂x2 = 1000 cos x2 - 0.25 + 1000 cos x2 - x1 - 0.25 ∂f6 ∂x3 = 3×10 -6 x32 ∂f6 ∂x4 = 2×10 -6 x42
So the arrays igfun and jgvar must contain:
k 1 2 3 4 5 6 7 8
igfun[k-1] 1 1 2 2 3 3 6 6
jgvar[k-1] 1 2 1 2 1 2 3 4
and usrfun must return in g[k-1]  the value of ∂fi ∂xj , where i=igfun[k-1]  and j=jgvar[k-1] .

Example program (C#): e04vhe.cs

Example program data: e04vhe.d

Example program results: e04vhe.r

Algorithmic Details

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

Constraints and Slack Variables

Problem (1) contains n  variables in x . Let m  be the number of components of fx  and ALx  combined. The upper and lower bounds on those terms define the general constraints of the problem. e04vh converts the general constraints to equalities by introducing a set of slack variables s = s1,s2,…,smT . 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 f0x  subject to ​ fx ALx - s = 0 ,  l ≤ x s ≤ u . (3)
The general constraints become the equalities fx - sN = 0  and ALx - sL = 0 , where sL  and sN  are the linear and nonlinear slacks.

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 . The constraints of the subproblem are formed from the linear constraints ALx - sL = 0  and the linearized constraint
fxk + f′ xk x-xk - sN = 0 , (4)
where f′ xk  denotes the Jacobian matrix, whose elements are the first derivatives of fx  evaluated at xk . The QP constraints therefore comprise the m  linear constraints
f′ xk x - sN = -fxk + f′ xk xk , ALx -sN - sL = 0 , (5)
where x  and s  are bounded above and below by u  and l  as before. If the m  by n  matrix A  and m -vector b  are defined as
A = f′ xk AL  and  b = -fxk + f′ xk xk 0 , (6)
then the QP subproblem can be written as
minimize x,s qx,xk = gkT x-xk + 12 x-xkT Hk x-xk  subject to ​ Ax - s = b ,  l ≤ x s ≤ u , (7)
where qx,xk  is a quadratic approximation to a modified Lagrangian function (see Gill et al. (2002)). 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 n1  rows and columns of the Hessian need to be approximated, where n1  is the number of nonlinear variables. This quantity is determined by the implicit values of the number of nonlinear objective and Jacobian variables determined after the constraints and variables are reordered.

Minor Iterations

Solving the QP subproblem is itself an iterative procedure. Here, the iterations of the QP solver e04nq form the minor iterations of the SQP method. e04nq 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 , (8)
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 g-ATπ , where g  is the gradient of the QP objective, and the reduced gradients for the slacks are the dual variables π . The QP subproblem is optimal if dj≥0  for all nonbasic variables at their lower bounds, dj≤0  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.

The Merit Function

Treatment of Constraint Infeasibilities

e04vh 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 ALx - v + w ≤ u ,  v≥0 ,  w≥0 , (11)
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 (11) has v≠0  or w≠0 . e04vh 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 functions can be safely evaluated.) e04vh proceeds to solve nonlinear problems as given, using search directions obtained from the sequence of QP subproblems (see (7)).
If a QP subproblem proves to be infeasible or unbounded (or if the dual variables π  for the nonlinear constraints become large), e04vh enters ‘elastic’ mode and thereafter solves the problem
minimize x,v,w f0x + γ eT v+w  subject to ​ l ≤ x fx - v + w ALx ≤ u ,  v≥0 ,  w≥0 , (12)
where γ  is a non-negative parameter (the elastic weight), and f0x + γ 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.
The initial value of γ  is controlled by the optional parameter Elastic Weight.

Description of Monitoring Information

e04vh produces monitoring information, statistical information and information about the solution. [The Final Output] contains details of the final output information sent to the unit specified by the optional parameter Print File. This section contains other details of output information.

Major Iteration Log

This section describes the output to unit Print File if Major Print Level > 0 . One line of information is output every k th major iteration, where k  is Print Frequency.
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 [Algorithmic Details]).
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 usrfun has 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 (see (13)), the maximum component of the scaled nonlinear constraint residual (see optional parameter Major Feasibility Tolerance). The solution is regarded as acceptably feasible if Feasible is less than the Major Feasibility Tolerance. 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 (see (14)), the maximum complementary gap (see optional parameter Major Optimality Tolerance). 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 (see (8)). This function will decrease at each iteration unless it was necessary to increase the penalty parameters (see [The Merit Function]). 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 (see [Basis Factorization Statistics]) is the number of subdiagonal elements in the columns of a lower triangular matrix and lenU (see [Basis Factorization Statistics]) is the number of diagonal and superdiagonal elements in the rows of an upper-triangular matrix.
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.
If the constraints are linear, refactorization is subject only to the Factorization Frequency, and L+U will tend to increase between factorizations.
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 follow the separator ‘_’ in the output
Label Description
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 conforming to the value of the optional parameter Violation Limit.
D you set status=-1  on exit from usrfun, indicating that the linesearch needed to be done with a smaller value of the step length α.
l the norm wise change in the variables was limited by the value of the optional parameter Major Step Limit. If this output occurs repeatedly during later iterations, it may be worthwhile increasing the value of the optional parameter Major Step Limit.
i if e04vh 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 (12) (see [Treatment of Constraint Infeasibilities]).
If e04vh is already in elastic mode, an i indicates that the minimizer of the elastic subproblem does not satisfy the linearized constraints. (In this case, a feasible point for the usual QP subproblem may or may not exist.)
M an extra evaluation of the problem functions was needed to define an acceptable positive-definite quasi-Newton update to the Lagrangian Hessian. This modification is only done when there are nonlinear constraints.
m this is the same as M except that it was also necessary to modify the update to include an augmented Lagrangian term.
n no positive-definite BFGS update could be found. The approximate Hessian is unchanged from the previous iteration.
R the approximate Hessian has been reset by discarding all but the diagonal elements. This reset will be forced periodically by the Hessian Frequency and Hessian Updates keywords. However, it may also be necessary to reset an ill-conditioned Hessian from time to time.
r the approximate Hessian was reset after ten consecutive major iterations in which no BFGS update could be made. The diagonals of the approximate Hessian are retained if at least one update has been done since the last reset. Otherwise, the approximate Hessian is reset to the identity matrix.
s a self-scaled BFGS update was performed. This update is used when the Hessian approximation is diagonal, and hence always follows a Hessian reset.
t the minor iterations were terminated because of the Minor Iterations Limit.
T the minor iterations were terminated because of the New Superbasics Limit.
u the QP subproblem was unbounded.
w a weak solution of the QP subproblem was found.
z the Superbasics Limit was reached.

Minor Iteration Log

If Minor Print Level > 0 , one line of information is output to the Print file every k th minor iteration, where k  is the specified Print Frequency. 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.
If Partial Price is in effect, variable jq is selected from App  or Ipp , the pp th segments of the constraint matrix A -I .
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 one 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). In rare cases, it may be necessary to increase the Pivot Tolerance to exclude very small elements of y  from consideration during the computation of Step.
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 [Major Iteration Log]. (The heading is not printed if the problem is linear.)

Crash Statistics

Basis Factorization Statistics

If Major Print Level ≥ 10 , the first seven items listed below are output to the Print file whenever the basis B  or the rectangular matrix BS = B S T  is factorized before solution of the next QP subproblem (see [Description of the Optional Parameters]).
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 . Stability is ensured as described under optional parameter LU Factor Tolerance.
If Minor Print Level ≥ 10 , the same items are printed during the QP solution whenever the current B  is factorized. In addition, if System Information Yes has been specified, the entries from Elems onwards are also printed.
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 Factorization Frequency.
2 The nonzeros in the updated factors have increased significantly.
7 Not enough storage to update factors.
10 Row residuals are too large (see the description of the optional parameter Check Frequency).
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. If LU Complete Pivoting or LU Rook Pivoting has been selected, this heading is changed to MerCP, respectively 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 e04vh 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 . This is the specified LU Factor Tolerance or a smaller value that is currently being used for greater stability.
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. If it is orders of magnitude larger, it may be advisable to reduce the LU Factor Tolerance to 5.0, 4.0, 3.0 or 2.0, say (but bigger than 1.0).
As long as Lmax is not large (say 5.0 or less), maxAmax,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.) Akmax is not printed if LU Partial Pivoting is selected.
Agrwth is the ratio Akmax/Amax . Values much larger than 100 (say) indicate instability. Agrwth is not printed if LU Partial Pivoting is selected.
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 of PUQ .
DUmin is the smallest diagonal 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).

The Summary File

See Also