naginterfaces.library.mip.sqp

naginterfaces.library.mip.sqp(ncnln, bl, bu, varcon, x, objfun, comm, a=None, d=None, confun=None, maxit=500, acc=0.0001, data=None, io_manager=None, spiked_sorder='C')

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

Note: this function uses optional algorithmic parameters, see also: optset().

For full information please refer to the NAG Library document for h02da

https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/h/h02daf.html

Parameters

ncnlnint

$n_N$, the number of constraints supplied by $c\left(x\right)$.

blfloat, array-like, shape $\left(n\right)$

$bl$ must contain the lower bounds, $l_i$, and $bu$ the upper bounds, $u_i$, for the variables; bounds on integer variables are rounded, bounds on binary variables need not be supplied.

bufloat, array-like, shape $\left(n\right)$

$bl$ must contain the lower bounds, $l_i$, and $bu$ the upper bounds, $u_i$, for the variables; bounds on integer variables are rounded, bounds on binary variables need not be supplied.

varconint, array-like, shape $(n+\text{nclin}+ncnln)$

$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 $n_L$ elements the nature of the general linear constraints (if any) and the next $n_N$ elements the general constraints (if any).

$varcon[j-1]=0$

A continuous variable.

$varcon[j-1]=1$

A binary variable.

$varcon[j-1]=2$

An integer variable.

$varcon[j-1]=3$

An equality constraint.

$varcon[j-1]=4$

An inequality constraint.

xfloat, array-like, shape $\left(n\right)$

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.

objfuncallable (objmip, objgrd) = objfun(mode, varcon, x, objgrd, nstate, data=None)

$objfun$ must calculate the objective function $f\left(x\right)$ and its gradient for a specified $n$-element vector $x$.

Parameters

modeint

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$.

varconint, ndarray, shape $\left(\text{lvarcon}\right)$

The array $varcon$ as supplied to sqp.

xfloat, ndarray, shape $\left(n\right)$

The vector of variables at which the objective function and/or all continuous elements of its gradient are to be evaluated.

objgrdfloat, ndarray, shape $\left(n\right)$

Continuous elements of $objgrd$ are set to the value of NaN.

nstateint

If $nstate=1$, sqp is calling $objfun$ for the first time. This argument setting allows you to save computation time if certain data must be read or calculated only once.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

objmipfloat

Must be set to the objective function value, $f$.

objgrdfloat, array-like, shape $\left(n\right)$

Must contain the gradient vector of the objective function if $mode=1$, with $objgrd[j-1]$ containing the partial derivative of $f$ with respect to continuous variable $x_j$; otherwise $objgrd$ is not referenced.

commdict, communication object

Communication structure.

This argument must have been initialized by a prior call to optset().

aNone or float, array-like, shape $(\text{nclin},:)$, optional

Note: the required extent for this argument in dimension 2 is determined as follows: if $\text{nclin}>0$: $n$; otherwise: $0$.

The $\text{i}$th row of $a$ must contain the coefficients of the $\text{i}$th general linear constraint, for $\text{i}=1,2,\dots ,n_l$. Any equality constraints must be specified first.

If $\text{nclin}=0$, the array $a$ is not referenced.

dNone or float, array-like, shape $\left(\text{nclin}\right)$, optional

$d_i$, the constant for the $i$th linear constraint.

If $\text{nclin}=0$, the array $d$ is not referenced.

confunNone or callable (c, cjac) = confun(mode, varcon, x, cjac, nstate, data=None), optional

Note: if this argument is None then a NAG-supplied facility will be used.

$confun$ must calculate the constraint functions supplied by $c\left(x\right)$ and their Jacobian at $x$. If all constraints are supplied by $A$ and $d$ (i.e., $ncnln=0$), $confun$ will never be called by sqp and $confun$ may be None. If $ncnln>0$, the first call to $confun$ will occur after the first call to $objfun$.

Parameters

modeint

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.

varconint, ndarray, shape $\left(\text{lvarcon}\right)$

The array $varcon$ as supplied to sqp.

xfloat, ndarray, shape $\left(n\right)$

The vector of variables at which the objective function and/or all continuous elements of its gradient are to be evaluated.

cjacfloat, ndarray, shape $(\text{ncnln},n)$

Note: the derivative of the $i$th constraint with respect to the $j$th variable, $\frac{\partial c_i}{\partial x_j}$, is stored in $cjac[i-1,j-1]$.

Continuous elements of $cjac$ are set to the value of NaN.

nstateint

If $nstate=1$, sqp is calling $confun$ for the first time. This argument setting allows you to save computation time if certain data must be read or calculated only once.

dataarbitrary, optional, modifiable in place

User-communication data for callback functions.

Returns

cfloat, array-like, shape $\left(\text{ncnln}\right)$

Must contain $ncnln$ constraint values, with the value of the $j$th constraint $c_j\left(x\right)$ in $c[j-1]$.

cjacfloat, array-like, shape $(\text{ncnln},n)$

The $i$th row of $cjac$ must contain elements of $\frac{\partial c_i}{\partial x_j}$ for each continuous variable $x_j$.

maxitint, optional

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.

accfloat, optional

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 machine.precision), the below suggested value is used.

dataarbitrary, optional

User-communication data for callback functions.

io_managerFileObjManager, optional

Manager for I/O in this routine.

spiked_sorderstr, optional

If $a$ is spiked (i.e., has unit extent in all but one dimension, or has size $1$), $spiked\_sorder$ selects the storage order to associate with it in the NAG Engine:

spiked_sorder = $\text{'C'}$

row-major storage will be used;

spiked_sorder = $\text{'F'}$

column-major storage will be used.

Two-dimensional arrays returned from callback functions in this routine must then use the same storage order.

Returns

axfloat, ndarray, shape $\left(\text{nclin}\right)$

The final residuals of the linear constraints $Ax-d$.

If $\text{nclin}=0$, $ax$ is not referenced.

xfloat, ndarray, shape $\left(n\right)$

The final estimate of the solution.

cfloat, ndarray, shape $\left(ncnln\right)$

If $ncnln>0$, $c[\text{j}-1]$ contains the value of the $\text{j}$th constraint function $c_{\text{j}}\left(x\right)$ at the final iterate, for $\text{j}=1,2,\dots ,ncnln$.

If $ncnln=0$, the array $c$ is not referenced.

cjacfloat, ndarray, shape $(ncnln,n)$

If $ncnln>0$, $cjac$ contains the Jacobian matrix of the constraint functions at the final iterate, i.e., $cjac[\text{i}-1,\text{j}-1]$ contains the partial derivative of the $\text{i}$th constraint function with respect to the $\text{j}$th variable, for $\text{j}=1,2,\dots ,n$, for $\text{i}=1,2,\dots ,ncnln$. (See the discussion of argument $cjac$ under $confun$.)

If $ncnln=0$, the array $cjac$ is not referenced.

objgrdfloat, ndarray, shape $\left(n\right)$

The objective function gradient at the solution.

objmipfloat

With no exception or warning is raised, $objmip$ contains the value of the objective function for the MINLP solution.

Other Parameters

‘Branch Bound Steps’int

Default $=500$

Maximum number of branch-and-bound steps for solving the mixed integer quadratic problems.

‘Branching Rule’str

Default $=\text{'Maximum'}$

Branching rule for branch and bound search.

$\text{‘Branching Rule'}=\text{'Maximum'}$

Maximum fractional branching.

$\text{‘Branching Rule'}=\text{'Minimum'}$

Minimum fractional branching.

‘Check Gradients’str

Default $=\text{'No'}$

Perform an internal check of supplied objective and constraint gradients. It is advisable to set $\text{‘Check Gradients'}=\text{'Yes'}$ during code development to avoid difficulties associated with incorrect user-supplied gradients.

‘Continuous Trust Radius’float

Default $=10.0$

Initial continuous trust region radius, ${\Delta }_0^c$; the initial trial step $d\in R^{n_c}$ for the SQP approximation must satisfy ${\parallel d\parallel }_{\infty }\le {\Delta }_0^c$.

‘Descent’float

Default $=0.05$

Initial descent parameter, $\delta$, larger values of $\delta$ allow penalty option $\sigma$ to increase faster which can lead to faster convergence.

‘Descent Factor’float

Default $=0.1$

Factor for decreasing the internal descent parameter, $\delta$, between iterations.

‘Feasible Steps’int

Default $=10$

Maximum number of feasible steps without improvements, where feasibility is measured by ${\parallel g\left(x\right)\parallel }_{\infty }\le \sqrt{acc}$.

‘Improved Bounds’str

Default $=\text{'No'}$

Calculate improved bounds in case of ‘Best of all’ node selection strategy.

‘Integer Trust Radius’float

Default $=10.0$

Initial integer trust region radius, ${\Delta }_0^i$; the initial trial step $e\in R^{n_i}$ for the SQP approximation must satisfy ${\parallel e\parallel }_{\infty }\le {\Delta }_0^i$.

‘Maximum Restarts’int

Default $=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.

‘Minor Print Level’int

Default $=0$

Print level of the subproblem solver. Active only if $\text{‘Print Level'}\ne 0$.

‘Modify Hessian’str

Default $=\text{'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 Selection’str

Default $=\text{'Depth First'}$

Node selection strategy for branch and bound.

$\text{‘Node Selection'}=\text{'Best of all'}$

Large tree search; this method is the slowest as it solves all subproblem QPs independently.

$\text{‘Node Selection'}=\text{'Best of two'}$

Uses warm starts and less memory.

$\text{‘Node Selection'}=\text{'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 Monotone’int

Default $=10$

Maximum number of successive iterations considered for the non-monotone trust region algorithm, allowing the penalty function to increase.

‘Objective Scale Bound’float

Default $=1.0$

When $\text{‘Scale Objective Function'}>0$ internally scale absolute function values greater than $1.0$ or ‘Objective Scale Bound’.

‘Penalty’float

Default $=1000.0$

Initial penalty option, $\sigma$.

‘Penalty Factor’float

Default $=10.0$

Factor for increasing penalty option $\sigma$ when the trust regions cannot be enlarged at a trial step.

‘Print Level’int

Default $=0$

Specifies the desired output level of printing.

$\text{‘Print Level'}=\text{'0'}$

No output.

$\text{‘Print Level'}=\text{'1'}$

Final convergence analysis.

$\text{‘Print Level'}=\text{'2'}$

One line of intermediate results per iteration.

$\text{‘Print Level'}=\text{'3'}$

Detailed information printed per iteration.

‘QP Accuracy’float

Default $=1.0e-10$

Termination tolerance of the relaxed quadratic program subproblems.

‘Scale Continuous Variables’str

Default $=\text{'Yes'}$

Internally scale continuous variables values.

‘Scale Objective Function’int

Default $=1$

Internally scale objective function values.

$\text{‘Scale Objective Function'}=\text{'0'}$

No scaling.

$\text{‘Scale Objective Function'}=\text{'1'}$

Scale absolute values greater than ‘Objective Scale Bound’.

‘Warm Starts’int

Default $=100$

Maximum number of warm starts within the mixed integer QP solver, see ‘Node Selection’.

Raises

NagValueError

(errno $1$)

On entry, $n=⟨value⟩$.

Constraint: $n>0$.

(errno $2$)

On entry, $\text{nclin}=⟨value⟩$.

Constraint: $\text{nclin}\ge 0$.

(errno $3$)

On entry, $ncnln=⟨value⟩$.

Constraint: $ncnln\ge 0$.

(errno $5$)

On entry, $bl[⟨value⟩]>bu[⟨value⟩]$.

Constraint: $bl[\text{i}-1]\le bu[\text{i}-1]$, for $\text{i}=1,2,\dots ,n$.

(errno $6$)

On entry, $varcon[⟨value⟩]=⟨value⟩$.

Constraint: $varcon[\text{i}-1]=0$, $1$ or $2$, for $\text{i}=1,2,\dots ,n$.

(errno $7$)

On entry, $varcon[⟨value⟩]=⟨value⟩$.

Constraint: $varcon[\text{i}-1]=3$ or $4$, for $\text{i}=n+1,\dots ,n+\text{nclin}+ncnln$.

(errno $8$)

The supplied $objfun$ returned a NaN value.

(errno $9$)

The supplied $confun$ returned a NaN value.

(errno $10$)

On entry, the option arrays $comm$[‘iopts’] and $comm$[‘opts’] have either not been initialized or been corrupted.

(errno $11$)

On entry, there are no binary or integer variables.

(errno $12$)

On entry, linear equality constraints do not precede linear inequality constraints.

(errno $13$)

On entry, nonlinear equality constraints do not precede nonlinear inequality constraints.

(errno $81$)

One or more objective gradients appear to be incorrect.

(errno $91$)

One or more constraint gradients appear to be incorrect.

(errno $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. Option $\text{‘Feasible Steps'}=⟨value⟩$.

(errno $1003$)

Penalty parameter tends to infinity in an underlying mixed-integer quadratic program; the problem may be infeasible. If $\sigma$ is relatively low value, try a higher one, for example ${10}^{20}$.

Option $\text{‘Penalty'}=⟨value⟩$.

(errno $1004$)

Termination at an infeasible iterate; if the problem is feasible, try a different starting value.

(errno $1005$)

Termination with zero integer trust region for integer variables; try a different starting value.

Option $\text{‘Integer Trust Radius'}=⟨value⟩$.

(errno $1008$)

The optimization failed due to numerical difficulties. Set option $\text{‘Print Level'}=3$ for more information.

Warns

NagAlgorithmicWarning

(errno $1001$)

On entry, $maxit=⟨value⟩$. Exceeded the maximum number of iterations.

NagCallbackTerminateWarning

(errno $i<0$)

The optimization halted because you set $mode$ negative in $objfun$ or $mode$ negative in $confun$, to $⟨value⟩$.

Notes

sqp 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:

$$\begin{array}{c}\begin{array}{ccc}{minimize}_{x\in \{R^{n_c},Z^{n_i}\}}& f\left(x\right)& \\ \text{subject to}& c_j\left(x\right)=0\text{,}& j=1,2,\dots ,m_e\\ & c_j\left(x\right)\ge 0\text{,}& j=m_e+1,m_e+2,\dots ,m\\ & l\le x_i\le u\text{,}& i=1,2,\dots ,n\end{array}\end{array}$$

with $n_c$ continuous variables and $n_i$ binary and integer variables in a total of $n$ variables; $m_e$ equality constraints in a total of $m$ constraint functions.

Partial derivatives of $f\left(x\right)$ and $c\left(x\right)$ are not required for the $n_i$ integer variables. Gradients with respect to integer variables are approximated by difference formulae.

No assumptions are made regarding $f\left(x\right)$ 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_{\sigma }\left(x\right)=f\left(x\right)+\sigma {\parallel g\left(x\right)\parallel }_{\infty }$$

where $\sigma$ is adapted by the algorithm and $g\left(x\right)$ is given by:

$$\begin{array}{c}\begin{array}{cccc}g\left(x\right)& =& c_j\left(x\right)\text{,}& j=1,2,\dots ,m_e\\ & =& min(c_j\left(x\right),0)\text{,}& j=m_e+1,m_e+2,\dots ,m\text{.}\end{array}\end{array}$$

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 The Quasi-Newton Update for opt.nlp1_solve) 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 $A$ and vector $d$ rather than the constraint functions $c\left(x\right)$ such that

$$Ax=d\text{or}Ax\ge d\text{.}$$

Partial derivatives with respect to $x$ of these constraint functions are not requested by sqp.

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

See also

naginterfaces.library.examples.mip.sqp_ex.main()