H02DAF (PDF version)
H Chapter Contents
H Chapter Introduction
NAG Library Manual

NAG Library Routine Document

H02DAF

Note:  before using this routine, please read the Users' Note for your implementation to check the interpretation of bold italicised terms and other implementation-dependent details.

 Contents

    1  Purpose
    7  Accuracy

1  Purpose

H02DAF solves general nonlinear programming problems with integer constraints on some of the variables.

2  Specification

SUBROUTINE H02DAF ( N, NCLIN, NCNLN, A, LDA, D, AX, BL, BU, VARCON, X, CONFUN, C, CJAC, OBJFUN, OBJGRD, MAXIT, ACC, OBJMIP, IOPTS, OPTS, IUSER, RUSER, IFAIL)
INTEGER  N, NCLIN, NCNLN, LDA, VARCON(N+NCLIN+NCNLN), MAXIT, IOPTS(*), IUSER(*), IFAIL
REAL (KIND=nag_wp)  A(LDA,*), D(NCLIN), AX(NCLIN), BL(N), BU(N), X(N), C(NCNLN), CJAC(NCNLN,N), OBJGRD(N), ACC, OBJMIP, OPTS(*), RUSER(*)
EXTERNAL  CONFUN, OBJFUN
Before calling H02DAF, H02ZKF must be called with OPTSTR set to ‘Initialize = h02daf’. Optional parameters may also be specified by calling H02ZKF before the call to H02DAF.

3  Description

H02DAF solves mixed integer nonlinear programming problems using a modified sequential quadratic programming method. The problem is assumed to be stated in the following general form:
minimize x Rnc,Zni fx subject to cjx=0, j=1,2,,me cjx0, j=me+1,me+2,,m lxiu, i=1,2,,n  
with nc continuous variables and ni binary and integer variables in a total of n variables; me equality constraints in a total of m constraint functions.
Partial derivatives of fx and cx are not required for the ni integer variables. Gradients with respect to integer variables are approximated by difference formulae.
No assumptions are made regarding fx except that it is twice continuously differentiable with respect to continuous elements of x. It is not assumed that integer variables are relaxable. In other words, problem functions are evaluated only at integer points.
The method seeks to minimize the exact penalty function:
Pσx = fx +σ gx  
where σ is adapted by the algorithm and gx is given by:
gx = cjx, j=1,2,,me = mincjx,0, j=me+1,me+2,,m.  
Successive quadratic approximations are applied under the assumption that integer variables have a smooth influence on the model functions, that is function values do not change drastically when incrementing or decrementing an integer value. In practice this requires integer variables to be ordinal not categorical. The algorithm is stabilised by a trust region method including Yuan's second order corrections, see Yuan and Sun (2006). The Hessian of the Lagrangian function is approximated by BFGS (see Section 11.4 in E04UCF/E04UCA) updates subject to the continuous and integer variables.
The mixed-integer quadratic programming subproblems of the SQP-trust region method are solved by a branch and cut method with continuous subproblem solutions obtained by the primal-dual method of Goldfarb and Idnani, see Powell (1983). Different strategies are available for selecting a branching variable: and a node from where branching, that is the generation of two new subproblems, begins:
This implementation is based on the algorithm MISQP as described in Exler et al. (2013).
Linear constraints may optionally be supplied by a matrix A and vector d rather than the constraint functions cx such that
Ax=d   or   ​ Axd .  
Partial derivatives with respect to x of these constraint functions are not requested by H02DAF.

4  References

Exler O, Lehmann T and Schittkowski K (2013) A comparative study of SQP-type algorithms for nonlinear and nonconvex mixed-integer optimization Mathematical Programming Computation 4 383–412
Mann A (1986) GAMS/MINOS: Three examples Department of Operations Research Technical Report Stanford University
Powell M J D (1983) On the quadratic programming algorithm of Goldfarb and Idnani Report DAMTP 1983/Na 19 University of Cambridge, Cambridge
Yuan Y-x and Sun W (2006) Optimization Theory and Methods Springer–Verlag

5  Parameters

1:     N – INTEGERInput
On entry: n, the total number of variables, nc+ni.
Constraint: N>0.
2:     NCLIN – INTEGERInput
On entry: nl, the number of general linear constraints defined by A and d.
Constraint: NCLIN0.
3:     NCNLN – INTEGERInput
On entry: nN, the number of constraints supplied by cx.
Constraint: NCNLN0.
4:     ALDA* – REAL (KIND=nag_wp) arrayInput
Note: the second dimension of the array A must be at least N if NCLIN>0.
On entry: the ith row of A must contain the coefficients of the ith general linear constraint, for i=1,2,,nl. Any equality constraints must be specified first.
If NCLIN=0, the array A is not referenced.
5:     LDA – INTEGERInput
On entry: the first dimension of the array A as declared in the (sub)program from which H02DAF is called.
Constraint: LDAmax1,NCLIN.
6:     DNCLIN – REAL (KIND=nag_wp) arrayInput
On entry: di, the constant for the ith linear constraint.
If NCLIN=0, the array D is not referenced.
7:     AXNCLIN – REAL (KIND=nag_wp) arrayOutput
On exit: the final values of the linear constraints Ax.
If NCLIN=0, AX is not referenced.
8:     BLN – REAL (KIND=nag_wp) arrayInput
9:     BUN – REAL (KIND=nag_wp) arrayInput
On entry: BL must contain the lower bounds, li, and BU the upper bounds, ui, for the variables; bounds on integer variables are rounded, bounds on binary variables need not be supplied.
Constraint: BLiBUi, for i=1,2,,N.
10:   VARCONN+NCLIN+NCNLN – INTEGER arrayInput
On entry: VARCON indicates the nature of each variable and constraint in the problem. The first n elements of the array must describe the nature of the variables, the next nL elements the nature of the general linear constraints (if any) and the next nN elements the general constraints (if any).
VARCONj=0
A continuous variable.
VARCONj=1
A binary variable.
VARCONj=2
An integer variable.
VARCONj=3
An equality constraint.
VARCONj=4
An inequality constraint.
Constraints:
  • VARCONj=0, 1 or 2, for j=1,2,,N;
  • VARCONj=3 or 4, for j=N+1,,N+NCLIN+NCNLN;
  • At least one variable must be either binary or integer;
  • Any equality constraints must precede any inequality constraints.
11:   XN – REAL (KIND=nag_wp) arrayInput/Output
On entry: an initial estimate of the solution, which need not be feasible. Values corresponding to integer variables are rounded; if an initial value less than half is supplied for a binary variable the value zero is used, otherwise the value one is used.
On exit: the final estimate of the solution.
12:   CONFUN – SUBROUTINE, supplied by the NAG Library or the user.External Procedure
CONFUN must calculate the constraint functions supplied by cx and their Jacobian at x. If all constraints are supplied by A and d (i.e., NCNLN=0), CONFUN will never be called by H02DAF and CONFUN may be the dummy routine H02DDM. (H02DDM is included in the NAG Library.) If NCNLN>0, the first call to CONFUN will occur after the first call to OBJFUN.
The specification of CONFUN is:
SUBROUTINE CONFUN ( MODE, NCNLN, N, VARCON, X, C, CJAC, NSTATE, IUSER, RUSER)
INTEGER  MODE, NCNLN, N, VARCON(*), NSTATE, IUSER(*)
REAL (KIND=nag_wp)  X(N), C(NCNLN), CJAC(NCNLN,N), RUSER(*)
1:     MODE – INTEGERInput/Output
On entry: indicates which values must be assigned during each call of OBJFUN. Only the following values need be assigned:
MODE=0
Elements of C containing continuous variables.
MODE=1
Elements of CJAC containing continuous variables.
On exit: may be set to a negative value if you wish to terminate the solution to the current problem, and in this case H02DAF will terminate with IFAIL set to MODE.
2:     NCNLN – INTEGERInput
On entry: the dimension of the array C and the first dimension of the array CJAC as declared in the (sub)program from which H02DAF is called. The number of constraints supplied by cx, nN.
3:     N – INTEGERInput
On entry: the second dimension of the array CJAC as declared in the (sub)program from which H02DAF is called. n, the total number of variables, nc+ni.
4:     VARCON* – INTEGER arrayInput
On entry: the array VARCON as supplied to H02DAF.
5:     XN – REAL (KIND=nag_wp) arrayInput
On entry: the vector of variables at which the objective function and/or all continuous elements of its gradient are to be evaluated.
6:     CNCNLN – REAL (KIND=nag_wp) arrayOutput
On exit: must contain NCNLN constraint values, with the value of the jth constraint cjx in Cj.
7:     CJACNCNLNN – REAL (KIND=nag_wp) arrayInput/Output
Note: the derivative of the ith constraint with respect to the jth variable, ci xj , is stored in CJACij.
On entry: continuous elements of CJAC are set to the value of NaN.
On exit: the ith row of CJAC must contain elements of ci xj  for each continuous variable xj.
8:     NSTATE – INTEGERInput
On entry: if NSTATE=1, H02DAF is calling CONFUN for the first time. This parameter setting allows you to save computation time if certain data must be read or calculated only once.
9:     IUSER* – INTEGER arrayUser Workspace
10:   RUSER* – REAL (KIND=nag_wp) arrayUser Workspace
CONFUN must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which H02DAF is called. Parameters denoted as Input must not be changed by this procedure.
13:   CNCNLN – REAL (KIND=nag_wp) arrayOutput
On exit: if NCNLN>0, Cj contains the value of the jth constraint function cjx at the final iterate, for j=1,2,,NCNLN.
If NCNLN=0, the array C is not referenced.
14:   CJACNCNLNN – REAL (KIND=nag_wp) arrayOutput
Note: the derivative of the ith constraint with respect to the jth variable, ci xj , is stored in CJACij.
On exit: if NCNLN>0, CJAC contains the Jacobian matrix of the constraint functions at the final iterate, i.e., CJACij contains the partial derivative of the ith constraint function with respect to the jth variable, for i=1,2,,NCNLN and j=1,2,,N. (See the discussion of parameter CJAC under CONFUN.)
If NCNLN=0, the array CJAC is not referenced.
15:   OBJFUN – SUBROUTINE, supplied by the user.External Procedure
OBJFUN must calculate the objective function fx and its gradient for a specified n-element vector x.
The specification of OBJFUN is:
SUBROUTINE OBJFUN ( MODE, N, VARCON, X, OBJMIP, OBJGRD, NSTATE, IUSER, RUSER)
INTEGER  MODE, N, VARCON(*), NSTATE, IUSER(*)
REAL (KIND=nag_wp)  X(N), OBJMIP, OBJGRD(N), RUSER(*)
1:     MODE – INTEGERInput/Output
On entry: indicates which values must be assigned during each call of OBJFUN. Only the following values need be assigned:
MODE=0
The objective function value, OBJMIP.
MODE=1
The continuous elements of OBJGRD.
On exit: may be set to a negative value if you wish to terminate the solution to the current problem, and in this case H02DAF will terminate with IFAIL set to MODE.
2:     N – INTEGERInput
On entry: n, the total number of variables, nc+ni.
3:     VARCON* – INTEGER arrayInput
On entry: the array VARCON as supplied to H02DAF.
4:     XN – REAL (KIND=nag_wp) arrayInput
On entry: the vector of variables at which the objective function and/or all continuous elements of its gradient are to be evaluated.
5:     OBJMIP – REAL (KIND=nag_wp)Output
On exit: must be set to the objective function value, f, if MODE=0; otherwise OBJMIP is not referenced.
6:     OBJGRDN – REAL (KIND=nag_wp) arrayInput/Output
On entry: continuous elements of OBJGRD are set to the value of NaN.
On exit: must contain the gradient vector of the objective function if MODE=1, with OBJGRDj containing the partial derivative of f with respect to continuous variable xj; otherwise OBJGRD is not referenced.
7:     NSTATE – INTEGERInput
On entry: if NSTATE=1, H02DAF is calling OBJFUN for the first time. This parameter setting allows you to save computation time if certain data must be read or calculated only once.
8:     IUSER* – INTEGER arrayUser Workspace
9:     RUSER* – REAL (KIND=nag_wp) arrayUser Workspace
OBJFUN is called with the parameters IUSER and RUSER as supplied to H02DAF. You are free to use the arrays IUSER and RUSER to supply information to OBJFUN as an alternative to using COMMON global variables.
OBJFUN must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which H02DAF is called. Parameters denoted as Input must not be changed by this procedure.
16:   OBJGRDN – REAL (KIND=nag_wp) arrayOutput
On exit: the objective function gradient at the solution.
17:   MAXIT – INTEGERInput
On entry: the maximum number of iterations within which to find a solution. If MAXIT is less than or equal to zero, the suggested value below is used.
Suggested value: MAXIT=500.
18:   ACC – REAL (KIND=nag_wp)Input
On entry: the requested accuracy for determining feasible points during iterations and for halting the method when the predicted improvement in objective function is less than ACC. If ACC is less than or equal to ε (ε being the machine precision as given by X02AJF), the below suggested value is used.
Suggested value: ACC=0.0001.
19:   OBJMIP – REAL (KIND=nag_wp)Output
On exit: with IFAIL=0, OBJMIP contains the value of the objective function for the MINLP solution.
20:   IOPTS* – INTEGER arrayCommunication Array
Note: the dimension of this array is dictated by the requirements of associated functions that must have been previously called. This array must be the same array passed as argument IOPTS in the previous call to H02ZKF.
21:   OPTS* – REAL (KIND=nag_wp) arrayCommunication Array
Note: the dimension of this array is dictated by the requirements of associated functions that must have been previously called. This array must be the same array passed as argument OPTS in the previous call to H02ZKF.
On entry: the real option array as returned by H02ZKF.
22:   IUSER* – INTEGER arrayUser Workspace
23:   RUSER* – REAL (KIND=nag_wp) arrayUser Workspace
IUSER and RUSER are not used by H02DAF, but are passed directly to CONFUN and OBJFUN and may be used to pass information to these routines as an alternative to using COMMON global variables.
24:   IFAIL – INTEGERInput/Output
On entry: IFAIL must be set to 0, -1​ or ​1. If you are unfamiliar with this parameter you should refer to Section 3.3 in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value -1​ or ​1 is recommended. If the output of error messages is undesirable, then the value 1 is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is 0. 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:
IFAIL=1
On entry, N=value.
Constraint: N>0.
IFAIL=2
On entry, NCLIN=value.
Constraint: NCLIN0.
IFAIL=3
On entry, NCNLN=value.
Constraint: NCNLN0.
IFAIL=4
On entry, LDA=value and NCLIN=value.
Constraint: LDANCLIN.
IFAIL=5
On entry, BLvalue>BUvalue.
Constraint: BLiBUi, for i=1,2,,N.
IFAIL=6
On entry, VARCONvalue=value.
Constraint: VARCONi=0, 1 or 2, for i=1,2,,N.
IFAIL=7
On entry, VARCONvalue=value.
Constraint: VARCONi=3 or 4, for i=N+1,,N+NCLIN+NCNLN.
IFAIL=8
The supplied OBJFUN returned a NaN value.
IFAIL=9
The supplied CONFUN returned a NaN value.
IFAIL=10
On entry, the optional parameter arrays IOPTS and OPTS have either not been initialized or been corrupted.
IFAIL=11
On entry, there are no binary or integer variables.
IFAIL=12
On entry, linear equality constraints do not precede linear inequality constraints.
IFAIL=13
On entry, nonlinear equality constraints do not precede nonlinear inequality constraints.
IFAIL=81
One or more objective gradients appear to be incorrect.
IFAIL=91
One or more constraint gradients appear to be incorrect.
IFAIL=1001
On entry, MAXIT=value. Exceeded the maximum number of iterations.
IFAIL=1002
More than the maximum number of feasible steps without improvement in the objective function. If the maximum number of feasible steps is small, say less than 5, try increasing it. Optional parameter Feasible Steps=value.
IFAIL=1003
Penalty parameter tends to infinity in an underlying mixed-integer quadratic program; the problem may be infeasible. If σ is relatively low value, try a higher one, for example 1020.
Optional parameter Penalty=value.
IFAIL=1004
Termination at an infeasible iterate; if the problem is feasible, try a different starting value.
IFAIL=1005
Termination with zero integer trust region for integer variables; try a different starting value.
Optional parameter Integer Trust Radius=value.
IFAIL=1008
The optimization failed due to numerical difficulties. Set optional parameter Print Level=3 for more information.
IFAIL<0
The optimization halted because you set MODE negative in OBJFUN or MODE negative in CONFUN, to value.
IFAIL=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 3.8 in the Essential Introduction for further information.
IFAIL=-399
Your licence key may have expired or may not have been installed correctly.
See Section 3.7 in the Essential Introduction for further information.
IFAIL=-999
Dynamic memory allocation failed.
See Section 3.6 in the Essential Introduction for further information.

7  Accuracy

Not applicable.

8  Parallelism and Performance

H02DAF is not threaded by NAG in any implementation.
H02DAF makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9  Further Comments

None.

10  Example

Select a portfolio of at most p assets from n available with expected return ρ, is fully invested and that minimizes
xTΣx subject to  rTx = ρ i=1 n xi = 1 xi yi i=1 n yi p xi 0 yi = 0 or 1  
where
This example is taken from Mann (1986) with
r = 89127 Σ = 43-10 3610 -11100 0000 p = 3 ρ = 10.  
Linear constraints are supplied through both A and d, and CONFUN.

10.1  Program Text

Program Text (h02dafe.f90)

10.2  Program Data

None.

10.3  Program Results

Program Results (h02dafe.r)

11  Optional Parameters

This section can be skipped if you wish to use the default values for all optional parameters, otherwise, the following is a list of the optional parameters available and a full description of each optional parameter is provided in Section 11.1.

11.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 in order to return single options to their default states.
Keywords and character values are case insensitive, however they must be separated by at least one space.
H02ZKF can be called to supply options, one call being necessary for each optional parameter. For example,
Call H02ZKF('Check Gradients = Yes', iopts, liopts, opts, lopts, ifail)
H02ZKF should be consulted for a full description of the method of supplying optional parameters.
For H02DAF the maximum length of the parameter CVALUE used by H02ZLF is 12.
Branch Bound StepsiDefault=500
Maximum number of branch-and-bound steps for solving the mixed integer quadratic problems.
Constraint: Branch Bound Steps>1.
Branching RuleaDefault=Maximum
Branching rule for branch and bound search.
Branching Rule=Maximum
Maximum fractional branching.
Branching Rule=Minimum
Minimum fractional branching.
Check GradientsaDefault=No
Perform an internal check of supplied objective and constraint gradients. It is advisable to set Check Gradients=Yes during code development to avoid difficulties associated with incorrect user-supplied gradients.
Continuous Trust RadiusrDefault=10.0
Initial continuous trust region radius, Δ0c; the initial trial step dRnc for the SQP approximation must satisfy d Δ0c.
Constraint: Continuous Trust Radius>0.0.
DescentrDefault=0.05
Initial descent parameter, δ, larger values of δ allow penalty parameter σ to increase faster which can lead to faster convergence.
Constraint: 0.0<Descent<1.0.
Descent FactorrDefault=0.1
Factor for decreasing the internal descent parameter, δ, between iterations.
Constraint: 0.0<Descent Factor<1.0.
Feasible StepsiDefault=10
Maximum number of feasible steps without improvements, where feasibility is measured by gxACC.
Constraint: Feasible Steps>1.
Improved BoundsaDefault=No
Calculate improved bounds in case of ‘Best of all’ node selection strategy.
Integer Trust RadiusrDefault=10.0
Initial integer trust region radius, Δ0i; the initial trial step eRni for the SQP approximation must satisfy eΔ0i.
Constraint: Integer Trust Radius1.0.
Maximum RestartsiDefault=2
Maximum number of restarts that allow the mixed integer SQP algorithm to return to a better solution. Setting a value higher than the default might lead to better results at the expense of function evaluations.
Constraint: 0<Maximum Restarts15.
Minor Print LeveliDefault=0
Print level of the subproblem solver. Active only if Print Level0.
Constraint: 0<Minor Print Level<4.
Modify HessianaDefault=Yes
Modify the Hessian approximation in an attempt to get more accurate search directions. Calculation time is increased when the number of integer variables is large.
Node SelectionaDefault=Depth First
Node selection strategy for branch and bound.
Node Selection=Best of all
Large tree search; this method is the slowest as it solves all subproblem QPs independently.
Node Selection=Best of two
Uses warm starts and less memory.
Node Selection=Depth first
Uses more warm starts. If warm starts are applied, they can speed up the solution of mixed integer subproblems significantly when solving almost identical QPs.
Non MonotoneiDefault=10
Maximum number of successive iterations considered for the non-monotone trust region algorithm, allowing the penalty function to increase.
Constraint: 0<Non Monotone<100.
Objective Scale BoundrDefault=1.0
When Scale Objective Function>0 internally scale absolute function values greater than 1.0 or Objective Scale Bound.
Constraint: Objective Scale Bound>0.0.
PenaltyrDefault=1000.0
Initial penalty parameter, σ.
Constraint: Penalty0.0.
Penalty FactorrDefault=10.0
Factor for increasing penalty parameter σ when the trust regions cannot be enlarged at a trial step.
Constraint: Penalty Factor>1.0.
Print LeveliDefault=0
Specifies the desired output level of printing.
Print Level=0
No output.
Print Level=1
Final convergence analysis.
Print Level=2
One line of intermediate results per iteration.
Print Level=3
Detailed information printed per iteration.
QP AccuracyrDefault=1.0E−10
Termination tolerance of the relaxed quadratic program subproblems.
Constraint: QP Accuracy>0.0.
Scale Continuous VariablesaDefault=Yes
Internally scale continuous variables values.
Scale Objective FunctioniDefault=1
Internally scale objective function values.
Scale Objective Function=0
No scaling.
Scale Objective Function=1
Scale absolute values greater than Objective Scale Bound.
Warm StartsiDefault=100
Maximum number of warm starts within the mixed integer QP solver, see Node Selection.
Constraint: Warm Starts0.

H02DAF (PDF version)
H Chapter Contents
H Chapter Introduction
NAG Library Manual

© The Numerical Algorithms Group Ltd, Oxford, UK. 2015