NAG Library Routine Document
H02DAF
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:
with
continuous variables and
binary and integer variables in a total of
variables;
equality constraints in a total of
constraint functions.
Partial derivatives of and are not required for the integer variables. Gradients with respect to integer variables are approximated by difference formulae.
No assumptions are made regarding except that it is twice continuously differentiable with respect to continuous elements of . 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:
where
is adapted by the algorithm and
is given by:
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:
- Maximal fractional branching. Select an integer variable from the relaxed subproblem solution with largest distance from next integer value
- Minimal fractional branching. Select an integer variable from the relaxed subproblem solution with smallest distance from next integer value
and a node from where branching, that is the generation of two new subproblems, begins:
- Best of two. The optimal objective function values of the two child nodes are compared and the node with a lower value is chosen
- Best of all. Select an integer variable from the relaxed subproblem solution with the smallest distance from the next integer value
- Depth first. Select a child node whenever possible.
This implementation is based on the algorithm MISQP as described in
Exler et al. (2013).
Linear constraints may optionally be supplied by a matrix
and vector
rather than the constraint functions
such that
Partial derivatives with respect to 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: – INTEGERInput
-
On entry: , the total number of variables, .
Constraint:
.
- 2: – INTEGERInput
-
On entry: , the number of general linear constraints defined by and .
Constraint:
.
- 3: – INTEGERInput
-
On entry: , the number of constraints supplied by .
Constraint:
.
- 4: – REAL (KIND=nag_wp) arrayInput
-
Note: the second dimension of the array
A
must be at least
if
.
On entry: the
th row of
A must contain the coefficients of the
th general linear constraint, for
. Any equality constraints must be specified first.
If
, the array
A is not referenced.
- 5: – INTEGERInput
-
On entry: the first dimension of the array
A as declared in the (sub)program from which H02DAF is called.
Constraint:
.
- 6: – REAL (KIND=nag_wp) arrayInput
-
On entry:
, the constant for the
th linear constraint.
If
, the array
D is not referenced.
- 7: – REAL (KIND=nag_wp) arrayOutput
-
On exit: the final values of the linear constraints
.
If
,
AX is not referenced.
- 8: – REAL (KIND=nag_wp) arrayInput
- 9: – REAL (KIND=nag_wp) arrayInput
-
On entry: must contain the lower bounds, , and the upper bounds, , for the variables; bounds on integer variables are rounded, bounds on binary variables need not be supplied.
Constraint:
, for .
- 10: – INTEGER arrayInput
-
On entry:
VARCON indicates the nature of each variable and constraint in the problem. The first
elements of the array must describe the nature of the variables, the next
elements the nature of the general linear constraints (if any) and the next
elements the general constraints (if any).
- A continuous variable.
- A binary variable.
- An integer variable.
- An equality constraint.
- An inequality constraint.
Constraints:
- , or , for ;
- or , for ;
- At least one variable must be either binary or integer;
- Any equality constraints must precede any inequality constraints.
- 11: – 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: – SUBROUTINE, supplied by the NAG Library or the user.External Procedure
-
CONFUN must calculate the constraint functions supplied by
and their Jacobian at
. If all constraints are supplied by
and
(i.e.,
),
CONFUN will never be called by H02DAF and
CONFUN may be the dummy routine H02DDM. (H02DDM is included in the NAG Library.) If
, the first call to
CONFUN will occur after the first call to
OBJFUN.
The specification of
CONFUN is:
INTEGER |
MODE, NCNLN, N, VARCON(*), NSTATE, IUSER(*) |
REAL (KIND=nag_wp) |
X(N), C(NCNLN), CJAC(NCNLN,N), RUSER(*) |
|
- 1: – INTEGERInput/Output
-
On entry: indicates which values must be assigned during each call of
OBJFUN. Only the following values need be assigned:
- Elements of C containing continuous variables.
- 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: – 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
,
.
- 3: – INTEGERInput
-
On entry: the second dimension of the array
CJAC as declared in the (sub)program from which H02DAF is called.
, the total number of variables,
.
- 4: – INTEGER arrayInput
-
On entry: the array
VARCON as supplied to H02DAF.
- 5: – 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: – REAL (KIND=nag_wp) arrayOutput
-
On exit: must contain
NCNLN constraint values, with the value of the
th constraint
in
.
- 7: – REAL (KIND=nag_wp) arrayInput/Output
-
Note: the derivative of the th constraint with respect to the th variable, , is stored in .
On entry: continuous elements of
CJAC are set to the value of NaN.
On exit: the
th row of
CJAC must contain elements of
for each continuous variable
.
- 8: – INTEGERInput
-
On entry: if
, 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: – INTEGER arrayUser Workspace
- 10: – 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: – REAL (KIND=nag_wp) arrayOutput
-
On exit: if
,
contains the value of the
th constraint function
at the final iterate, for
.
If
, the array
C is not referenced.
- 14: – REAL (KIND=nag_wp) arrayOutput
-
Note: the derivative of the th constraint with respect to the th variable, , is stored in .
On exit: if
,
CJAC contains the Jacobian matrix of the constraint functions at the final iterate, i.e.,
contains the partial derivative of the
th constraint function with respect to the
th variable, for
and
. (See the discussion of parameter
CJAC under
CONFUN.)
If
, the array
CJAC is not referenced.
- 15: – SUBROUTINE, supplied by the user.External Procedure
-
OBJFUN must calculate the objective function
and its gradient for a specified
-element vector
.
The specification of
OBJFUN is:
INTEGER |
MODE, N, VARCON(*), NSTATE, IUSER(*) |
REAL (KIND=nag_wp) |
X(N), OBJMIP, OBJGRD(N), RUSER(*) |
|
- 1: – INTEGERInput/Output
-
On entry: indicates which values must be assigned during each call of
OBJFUN. Only the following values need be assigned:
- The objective function value, OBJMIP.
- 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: – INTEGERInput
-
On entry: , the total number of variables, .
- 3: – INTEGER arrayInput
-
On entry: the array
VARCON as supplied to H02DAF.
- 4: – 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: – REAL (KIND=nag_wp)Output
-
On exit: must be set to the objective function value,
, if
; otherwise
OBJMIP is not referenced.
- 6: – 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
, with
containing the partial derivative of
with respect to continuous variable
; otherwise
OBJGRD is not referenced.
- 7: – INTEGERInput
-
On entry: if
, 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: – INTEGER arrayUser Workspace
- 9: – 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: – REAL (KIND=nag_wp) arrayOutput
-
On exit: the objective function gradient at the solution.
- 17: – 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:
.
- 18: – 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:
.
- 19: – REAL (KIND=nag_wp)Output
-
On exit: with
,
OBJMIP contains the value of the objective function for the MINLP solution.
- 20: – 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: – 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: – INTEGER arrayUser Workspace
- 23: – 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: – INTEGERInput/Output
-
On entry:
IFAIL must be set to
,
. 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
is recommended. If the output of error messages is undesirable, then the value
is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is
.
When the value is used it is essential to test the value of IFAIL on exit.
On exit:
unless the routine detects an error or a warning has been flagged (see
Section 6).
6 Error Indicators and Warnings
If on entry
or
, explanatory error messages are output on the current error message unit (as defined by
X04AAF).
Errors or warnings detected by the routine:
-
On entry, .
Constraint: .
-
On entry, .
Constraint: .
-
On entry, .
Constraint: .
-
On entry, and .
Constraint: .
-
On entry, .
Constraint: , for .
-
On entry, .
Constraint: , or , for .
-
On entry, .
Constraint: or , for .
-
The supplied
OBJFUN returned a NaN value.
-
The supplied
CONFUN returned a NaN value.
-
On entry, the optional parameter arrays
IOPTS and
OPTS have either not been initialized or been corrupted.
-
On entry, there are no binary or integer variables.
-
On entry, linear equality constraints do not precede linear inequality constraints.
-
On entry, nonlinear equality constraints do not precede nonlinear inequality constraints.
-
One or more objective gradients appear to be incorrect.
-
One or more constraint gradients appear to be incorrect.
-
On entry, . Exceeded the maximum number of iterations.
-
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 , try increasing it. Optional parameter .
-
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 .
Optional parameter .
-
Termination at an infeasible iterate; if the problem is feasible, try a different starting value.
-
Termination with zero integer trust region for integer variables; try a different starting value.
Optional parameter .
-
The optimization failed due to numerical difficulties. Set optional parameter for more information.
-
The optimization halted because you set
MODE negative in
OBJFUN or
MODE negative in
CONFUN, to
.
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 3.8 in the Essential Introduction for further information.
Your licence key may have expired or may not have been installed correctly.
See
Section 3.7 in the Essential Introduction for further information.
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.
None.
10 Example
Select a portfolio of at most
assets from
available with expected return
, is fully invested and that minimizes
where
- is a vector of proportions of selected assets
- is an indicator variable that describes if an asset is in or out
- is a vector of mean returns
- is the covariance matrix of returns.
This example is taken from
Mann (1986) with
Linear constraints are supplied through both
and
, 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:
- the keywords;
- a parameter value,
where the letters , and denote options that take character, integer and real values respectively.
All options accept the value 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
.
Branch Bound Steps | | Default |
Maximum number of branch-and-bound steps for solving the mixed integer quadratic problems.
Constraint:
.
Branching Rule | | Default |
Branching rule for branch and bound search.
- Maximum fractional branching.
- Minimum fractional branching.
Check Gradients | | Default |
Perform an internal check of supplied objective and constraint gradients. It is advisable to set during code development to avoid difficulties associated with incorrect user-supplied gradients.
Continuous Trust Radius | | Default |
Initial continuous trust region radius, ; the initial trial step for the SQP approximation must satisfy .
Constraint:
.
Initial descent parameter, , larger values of allow penalty parameter to increase faster which can lead to faster convergence.
Constraint:
.
Descent Factor | | Default |
Factor for decreasing the internal descent parameter, , between iterations.
Constraint:
.
Feasible Steps | | Default |
Maximum number of feasible steps without improvements, where feasibility is measured by .
Constraint:
.
Improved Bounds | | Default |
Calculate improved bounds in case of ‘Best of all’ node selection strategy.
Integer Trust Radius | | Default |
Initial integer trust region radius, ; the initial trial step for the SQP approximation must satisfy .
Constraint:
.
Maximum Restarts | | Default |
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:
.
Minor Print Level | | Default |
Print level of the subproblem solver. Active only if .
Constraint:
.
Modify Hessian | | Default |
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 Selection | | Default |
Node selection strategy for branch and bound.
- Large tree search; this method is the slowest as it solves all subproblem QPs independently.
- Uses warm starts and less memory.
- 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.
Maximum number of successive iterations considered for the non-monotone trust region algorithm, allowing the penalty function to increase.
Constraint:
.
Objective Scale Bound | | Default |
When
internally scale absolute function values greater than
or
Objective Scale Bound.
Constraint:
.
Initial penalty parameter, .
Constraint:
.
Penalty Factor | | Default |
Factor for increasing penalty parameter when the trust regions cannot be enlarged at a trial step.
Constraint:
.
Specifies the desired output level of printing.
- No output.
- Final convergence analysis.
- One line of intermediate results per iteration.
- Detailed information printed per iteration.
QP Accuracy | | Default |
Termination tolerance of the relaxed quadratic program subproblems.
Constraint:
.
Scale Continuous Variables | | Default |
Internally scale continuous variables values.
Scale Objective Function | | Default |
Internally scale objective function values.
- No scaling.
- Scale absolute values greater than Objective Scale Bound.
Maximum number of warm starts within the mixed integer QP solver, see
Node Selection.
Constraint:
.