slsqp – slsqp

public subroutine slsqp(m, meq, la, n, x, xl, xu, f, c, g, a, acc, iter, mode, w, l_w, sdat, ldat, alphamin, alphamax, tolf, toldf, toldx, max_iter_ls, nnls_mode, infinite_bound)

slsqp: sequential least squares programming to solve general nonlinear optimization problems

a nonlinear programming method with quadratic programming subproblems this subroutine solves the general nonlinear programming problem:

minimize

$f(x)$

subject to

$c_j (x) = 0$ , $j = 1,...,meq$
$c_j (x) \ge 0$ , $j = meq+1,...,m$
$xl_i \le x_i \le xu_i$ , $i = 1,...,n$

the algorithm implements the method of Han and Powell with BFGS-update of the b-matrix and L1-test function within the steplength algorithm.

Reference

Dieter Kraft: "A software package for sequential quadratic programming", DFVLR-FB 88-28, 1988

History

implemented by: Dieter Kraft, DFVLR oberpfaffenhofen
date: april - october, 1981.
December, 31-st, 1984.
March , 21-st, 1987, revised to fortran 77
March , 20-th, 1989, revised to ms-fortran
April , 14-th, 1989, hesse in-line coded
February, 28-th, 1991, fortran/2 version 1.04 accepts statement functions
March , 1-st, 1991, tested with salford ftn77/386 compiler vers 2.40 in protected mode
January , 2016, Refactored into modern Fortran by Jacob Williams

License

Note

f, c, g, a must all be set by the user before each call.

Arguments

Type	Intent	Attributes		Name
integer,	intent(in)		::	m	is the total number of constraints, $m \ge 0$
integer,	intent(in)		::	meq	is the number of equality constraints, $m_{eq} \ge 0$
integer,	intent(in)		::	la	see `a`, $la \ge \max(m,1)$
integer,	intent(in)		::	n	is the number of variables, $n \ge 1$
real(kind=wp),	intent(inout),	dimension(n)	::	x	`x()` stores the current iterate of the `n` vector `x` on entry `x()` must be initialized. on exit `x()` stores the solution vector `x` if `mode = 0`.
real(kind=wp),	intent(in),	dimension(n)	::	xl	`xl()` stores an n vector of lower bounds `xl` to `x`.
real(kind=wp),	intent(in),	dimension(n)	::	xu	`xu()` stores an n vector of upper bounds `xu` to `x`.
real(kind=wp),	intent(in)		::	f	is the value of the objective function.
real(kind=wp),	intent(in),	dimension(la)	::	c	`c()` stores the `m` vector `c` of constraints, equality constraints (if any) first. dimension of `c` must be greater or equal `la`, which must be greater or equal `max(1,m)`.
real(kind=wp),	intent(in),	dimension(n+1)	::	g	`g()` stores the `n` vector `g` of partials of the objective function; dimension of `g` must be greater or equal `n+1`.
real(kind=wp),	intent(in),	dimension(la,n+1)	::	a	the `la` by `n + 1` array `a()` stores the `m` by `n` matrix `a` of constraint normals. `a()` has first dimensioning parameter `la`, which must be greater or equal `max(1,m)`.
real(kind=wp),	intent(inout)		::	acc	`abs(acc)` controls the final accuracy. if `acc` < zero an exact linesearch is performed, otherwise an armijo-type linesearch is used.
integer,	intent(inout)		::	iter	prescribes the maximum number of iterations. on exit `iter` indicates the number of iterations.
integer,	intent(inout)		::	mode	mode controls calculation: reverse communication is used in the sense that the program is initialized by `mode = 0`; then it is to be called repeatedly by the user until a return with `mode /= abs(1)` takes place. if `mode = -1` gradients have to be calculated, while with `mode = 1` functions have to be calculated. mode must not be changed between subsequent calls of slsqp. evaluation modes: -1 : gradient evaluation, (`g` & `a`) ** 0 : on entry: initialization, (`f`, `g`, `c`, `a`), on exit: required accuracy for solution obtained 1 : function evaluation, (`f` & `c`) failure modes: 2 : number of equality constraints larger than `n` 3 : more than `3n` iterations in lsq subproblem * 4 : inequality constraints incompatible 5 : singular matrix `e` in lsq subproblem 6 : singular matrix `c` in lsq subproblem 7 : rank-deficient equality constraint subproblem hfti 8 : positive directional derivative for linesearch 9 : more than `iter` iterations in sqp >=10 **: working space `w` too small, `w` should be enlarged to `l_w=mode/1000`,
real(kind=wp),	intent(inout),	dimension(l_w)	::	w	`w()` is a one dimensional working space. the first `m+n+nn1/2` elements of `w` must not be changed between subsequent calls of slsqp. on return `w(1) ... w(m)` contain the multipliers associated with the general constraints, while `w(m+1) ... w(m+n(n+1)/2)` store the cholesky factor `ldl(t)` of the approximate hessian of the lagrangian columnwise dense as lower triangular unit matrix `l` with `d` in its 'diagonal' and `w(m+n(n+1)/2+n+2 ... w(m+n(n+1)/2+n+2+m+2n)` contain the multipliers associated with all constraints of the quadratic program finding the search direction to the solution `x`
integer,	intent(in)		::	l_w	the length of `w`, which should be at least: `(3n1+m)(n1+1)` for lsq `+(n1-meq+1)(mineq+2) + 2mineq` for lsi `+(n1+mineq)(n1-meq) + 2meq + n1` for lsei `+ n1n/2 + 2m + 3n + 3n1 + 1` for slsqpb with `mineq = m - meq + 2*n1` & `n1 = n+1`
type(slsqpb_data),	intent(inout)		::	sdat	data for slsqpb.
type(linmin_data),	intent(inout)		::	ldat	data for linmin.
real(kind=wp),	intent(in)		::	alphamin	min $\alpha$ for line search $0 < \alpha_{min} < \alpha_{max} \le 1$
real(kind=wp),	intent(in)		::	alphamax	max $\alpha$ for line search $0 < \alpha_{min} < \alpha_{max} \le 1$
real(kind=wp),	intent(in)		::	tolf	stopping criterion if $\|f\| < tolf$ then stop.
real(kind=wp),	intent(in)		::	toldf	stopping criterion if $\|f_{n+1} - f_n\| < toldf$ then stop.
real(kind=wp),	intent(in)		::	toldx	stopping criterion if $\|\|x_{n+1} - x_n\|\| < toldx$ then stop.
integer,	intent(in)		::	max_iter_ls	maximum number of iterations in the nnls problem
integer,	intent(in)		::	nnls_mode	which NNLS method to use
real(kind=wp),	intent(in)		::	infinite_bound	"infinity" for the upper and lower bounds. if `xl<=-infinite_bound` or `xu>=infinite_bound` then these bounds are considered nonexistant. If `infinite_bound=0` then `huge()` is used for this.

Called by

Help

Source Code

    subroutine slsqp(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,w,l_w, &
                     sdat,ldat,alphamin,alphamax,tolf,toldf,toldx,&
                     max_iter_ls,nnls_mode,infinite_bound)

    implicit none

    integer,intent(in) :: m                     !! is the total number of constraints, \( m \ge 0 \)
    integer,intent(in) :: meq                   !! is the number of equality constraints, \( m_{eq} \ge 0 \)
    integer,intent(in) :: la                    !! see `a`, \( la \ge \max(m,1) \)
    integer,intent(in) :: n                     !! is the number of variables, \( n \ge 1 \)
    real(wp),dimension(n),intent(inout) :: x    !! `x()` stores the current iterate of the `n` vector `x`
                                                !! on entry `x()` must be initialized. on exit `x()`
                                                !! stores the solution vector `x` if `mode = 0`.
    real(wp),dimension(n),intent(in) :: xl      !! `xl()` stores an n vector of lower bounds `xl` to `x`.
    real(wp),dimension(n),intent(in) :: xu      !! `xu()` stores an n vector of upper bounds `xu` to `x`.
    real(wp),intent(in) :: f                    !! is the value of the objective function.
    real(wp),dimension(la),intent(in) :: c      !! `c()` stores the `m` vector `c` of constraints,
                                                !! equality constraints (if any) first.
                                                !! dimension of `c` must be greater or equal `la`,
                                                !! which must be greater or equal `max(1,m)`.
    real(wp),dimension(n+1),intent(in) :: g     !! `g()` stores the `n` vector `g` of partials of the
                                                !! objective function; dimension of `g` must be
                                                !! greater or equal `n+1`.
    real(wp),dimension(la,n+1),intent(in) ::  a !! the `la` by `n + 1` array `a()` stores
                                                !! the `m` by `n` matrix `a` of constraint normals.
                                                !! `a()` has first dimensioning parameter `la`,
                                                !! which must be greater or equal `max(1,m)`.
    real(wp),intent(inout) :: acc   !! `abs(acc)` controls the final accuracy.
                                    !! if `acc` < zero an exact linesearch is performed,
                                    !! otherwise an armijo-type linesearch is used.
    integer,intent(inout) :: iter   !! prescribes the maximum number of iterations.
                                    !! on exit `iter` indicates the number of iterations.
    integer,intent(inout) :: mode   !! mode controls calculation:
                                    !!
                                    !! reverse communication is used in the sense that
                                    !! the program is initialized by `mode = 0`; then it is
                                    !! to be called repeatedly by the user until a return
                                    !! with `mode /= abs(1)` takes place.
                                    !! if `mode = -1` gradients have to be calculated,
                                    !! while with `mode = 1` functions have to be calculated.
                                    !! mode must not be changed between subsequent calls of [[slsqp]].
                                    !!
                                    !! **evaluation modes**:
                                    !!
                                    !! * ** -1 **: gradient evaluation, (`g` & `a`)
                                    !! * **  0 **: *on entry*: initialization, (`f`, `g`, `c`, `a`),
                                    !!   *on exit*: required accuracy for solution obtained
                                    !! * **  1 **: function evaluation, (`f` & `c`)
                                    !!
                                    !! **failure modes**:
                                    !!
                                    !! * ** 2 **: number of equality constraints larger than `n`
                                    !! * ** 3 **: more than `3*n` iterations in [[lsq]] subproblem
                                    !! * ** 4 **: inequality constraints incompatible
                                    !! * ** 5 **: singular matrix `e` in [[lsq]] subproblem
                                    !! * ** 6 **: singular matrix `c` in [[lsq]] subproblem
                                    !! * ** 7 **: rank-deficient equality constraint subproblem [[hfti]]
                                    !! * ** 8 **: positive directional derivative for linesearch
                                    !! * ** 9 **: more than `iter` iterations in sqp
                                    !! * ** >=10 **: working space `w` too small,
                                    !!   `w` should be enlarged to `l_w=mode/1000`,
    integer,intent(in) :: l_w       !! the length of `w`, which should be at least:
                                    !!
                                    !! * `(3*n1+m)*(n1+1)`                     **for lsq**
                                    !! * `+(n1-meq+1)*(mineq+2) + 2*mineq`     **for lsi**
                                    !! * `+(n1+mineq)*(n1-meq) + 2*meq + n1`   **for lsei**
                                    !! * `+ n1*n/2 + 2*m + 3*n + 3*n1 + 1`     **for slsqpb**
                                    !!
                                    !! with `mineq = m - meq + 2*n1` & `n1 = n+1`
    real(wp),dimension(l_w),intent(inout) :: w  !! `w()` is a one dimensional working space.
                                                !! the first `m+n+n*n1/2` elements of `w` must not be
                                                !! changed between subsequent calls of [[slsqp]].
                                                !! on return `w(1) ... w(m)` contain the multipliers
                                                !! associated with the general constraints, while
                                                !! `w(m+1) ... w(m+n(n+1)/2)` store the cholesky factor
                                                !! `l*d*l(t)` of the approximate hessian of the
                                                !! lagrangian columnwise dense as lower triangular
                                                !! unit matrix `l` with `d` in its 'diagonal' and
                                                !! `w(m+n(n+1)/2+n+2 ... w(m+n(n+1)/2+n+2+m+2n)`
                                                !! contain the multipliers associated with all
                                                !! constraints of the quadratic program finding
                                                !! the search direction to the solution `x*`
    type(slsqpb_data),intent(inout) :: sdat  !! data for [[slsqpb]].
    type(linmin_data),intent(inout) :: ldat  !! data for [[linmin]].
    real(wp),intent(in) :: alphamin  !! min \( \alpha \) for line search
                                     !! \( 0 < \alpha_{min} < \alpha_{max} \le 1 \)
    real(wp),intent(in) :: alphamax  !! max \( \alpha \) for line search
                                     !! \( 0 < \alpha_{min} < \alpha_{max} \le 1 \)
    real(wp),intent(in) :: tolf      !! stopping criterion if \( |f| < tolf \) then stop.
    real(wp),intent(in) :: toldf     !! stopping criterion if \( |f_{n+1} - f_n| < toldf \) then stop.
    real(wp),intent(in) :: toldx     !! stopping criterion if \( ||x_{n+1} - x_n|| < toldx \) then stop.
    integer,intent(in)  :: max_iter_ls !! maximum number of iterations in the [[nnls]] problem
    integer,intent(in)  :: nnls_mode !! which NNLS method to use
    real(wp),intent(in) :: infinite_bound !! "infinity" for the upper and lower bounds.
                                          !! if `xl<=-infinite_bound` or `xu>=infinite_bound`
                                          !! then these bounds are considered nonexistant.
                                          !! If `infinite_bound=0` then `huge()` is used for this.

    integer :: il , im , ir , is , iu , iv , iw , ix , mineq, n1
    real(wp) :: infBnd !! local copy of `infinite_bound`

    if (infinite_bound==zero) then
        infBnd = huge(one)
    else
        infBnd = abs(infinite_bound)
    end if

    ! check length of working arrays
    n1 = n + 1
    mineq = m - meq + n1 + n1
    il = (3*n1+m)*(n1+1) + (n1-meq+1)*(mineq+2) + 2*mineq + (n1+mineq)&
         *(n1-meq) + 2*meq + n1*n/2 + 2*m + 3*n + 4*n1 + 1
    im = max(mineq,n1-meq)
    if ( l_w<il ) then
        mode = 1000*max(10,il)
        mode = mode + max(10,im)
        iter = 0
        return
    end if
    if (meq>n) then
        ! note: calling lsq when meq>n is corrupting the
        ! memory in some way, so just catch this here.
        mode = 2
        iter = 0
        return
    end if

    ! prepare data for calling sqpbdy  -  initial addresses in w

    im = 1
    il = im + max(1,m)
    il = im + la
    ix = il + n1*n/2 + 1
    ir = ix + n
    is = ir + n + n + max(1,m)
    is = ir + n + n + la
    iu = is + n1
    iv = iu + n1
    iw = iv + n1

    sdat%n1 = n1

    call slsqpb(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,&
                w(ir),w(il),w(ix),w(im),w(is),w(iu),w(iv),w(iw),&
                sdat%t,sdat%f0,sdat%h1,sdat%h2,sdat%h3,sdat%h4,&
                sdat%n1,sdat%n2,sdat%n3,sdat%t0,sdat%gs,sdat%tol,sdat%line,&
                sdat%alpha,sdat%iexact,sdat%incons,sdat%ireset,sdat%itermx,&
                ldat,alphamin,alphamax,tolf,toldf,toldx,&
                max_iter_ls,nnls_mode,infBnd)

    end subroutine slsqp

slsqp Subroutine

Contents

Source Code

public subroutine slsqp(m, meq, la, n, x, xl, xu, f, c, g, a, acc, iter, mode, w, l_w, sdat, ldat, alphamin, alphamax, tolf, toldf, toldx, max_iter_ls, nnls_mode, infinite_bound)

Reference

History

License

Arguments

Calls

Called by

Source Code