This subroutine seeks the least value of a function of many variables, subject to general linear inequality constraints, by a trust region method that forms quadratic models by interpolation.
LINCOA solves the following optimization problem:
Minimize F(X(1),X(2),...X(N)) subject to:
A * X <= B
Usually there is much freedom in each new model after satisfying the interpolation conditions, which is taken up by minimizing the Frobenius norm of the change to the second derivative matrix of the model. One new function value is calculated on each iteration, usually at a point where the current model predicts a reduction in the least value so far of the objective function subject to the linear constraints. Alternatively, a new vector of variables may be chosen to replace an interpolation point that may be too far away for reliability, and then the new point does not have to satisfy the linear constraints.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | n |
the number of variables. must be at least 2. |
||
| integer, | intent(in) | :: | npt |
the number of interpolation conditions, which is required to be in the interval [N+2,(N+1)(N+2)/2]. Typical choices of the author are NPT=N+6 and NPT=2*N+1. Larger values tend to be highly inefficent when the number of variables is substantial, due to the amount of work and extra difficulty of adjusting more points. |
||
| integer, | intent(in) | :: | m |
the number of linear inequality constraints. |
||
| real(kind=wp), | intent(in), | dimension(ia,*) | :: | a |
a matrix whose columns are the constraint gradients, which are required to be nonzero. |
|
| integer, | intent(in) | :: | ia |
the first dimension of the array A, which must be at least N. |
||
| real(kind=wp), | intent(in), | dimension(*) | :: | b |
the vector of right hand sides of the constraints, the J-th constraint being that the scalar product of A(.,J) with X(.) is at most B(J). The initial vector X(.) is made feasible by increasing the value of B(J) if necessary. |
|
| real(kind=wp), | intent(inout), | dimension(*) | :: | x |
the vector of variables. Initial values of X(1),X(2),...,X(N) must be supplied. If they do not satisfy the constraints, then B is increased as mentioned above. X contains on return the variables that have given the least calculated F subject to the constraints. |
|
| real(kind=wp), | intent(in) | :: | rhobeg |
RHOBEG and RHOEND must be set to the initial and final values of a trust region radius, so both must be positive with RHOEND<=RHOBEG. Typically, RHOBEG should be about one tenth of the greatest expected change to a variable, and RHOEND should indicate the accuracy that is required in the final values of the variables. |
||
| real(kind=wp), | intent(in) | :: | rhoend |
RHOBEG and RHOEND must be set to the initial and final values of a trust region radius, so both must be positive with RHOEND<=RHOBEG. Typically, RHOBEG should be about one tenth of the greatest expected change to a variable, and RHOEND should indicate the accuracy that is required in the final values of the variables. |
||
| integer, | intent(in) | :: | iprint |
The value of IPRINT should be set to 0, 1, 2 or 3, which controls the amount of printing. Specifically, there is no output if IPRINT=0 and there is output only at the return if IPRINT=1. Otherwise, the best feasible vector of variables so far and the corresponding value of the objective function are printed whenever RHO is reduced, where RHO is the current lower bound on the trust region radius. Further, each new value of F with its variables are output if IPRINT=3. |
||
| integer, | intent(in) | :: | maxfun |
an upper bound on the number of calls of CALFUN, its value being at least NPT+1. |
||
| procedure(func) | :: | calfun |
It must set F to the value of the objective function for the variables X(1), X(2),...,X(N). The value of the argument F is positive when CALFUN is called if and only if the current X satisfies the constraints |
subroutine lincoa (n, npt, m, a, ia, b, x, rhobeg, rhoend, iprint, maxfun, calfun) implicit none integer,intent(in) :: n !! the number of variables. must be at least 2. integer,intent(in) :: npt !! the number of interpolation conditions, which is !! required to be in the interval [N+2,(N+1)(N+2)/2]. Typical choices !! of the author are NPT=N+6 and NPT=2*N+1. Larger values tend to be !! highly inefficent when the number of variables is substantial, due !! to the amount of work and extra difficulty of adjusting more points. integer,intent(in) :: m !! the number of linear inequality constraints. integer,intent(in) :: ia !! the first dimension of the array A, which must be at least N. real(wp),dimension(ia,*),intent(in) :: a !! a matrix whose columns are the constraint gradients, which are !! required to be nonzero. real(wp),dimension(*),intent(in) :: b !! the vector of right hand sides of the constraints, the J-th !! constraint being that the scalar product of A(.,J) with X(.) is at !! most B(J). The initial vector X(.) is made feasible by increasing !! the value of B(J) if necessary. real(wp),dimension(*),intent(inout) :: x !! the vector of variables. Initial values of X(1),X(2),...,X(N) !! must be supplied. If they do not satisfy the constraints, then B !! is increased as mentioned above. X contains on return the variables !! that have given the least calculated F subject to the constraints. real(wp),intent(in) :: rhobeg !! RHOBEG and RHOEND must be set to the initial and final values of a !! trust region radius, so both must be positive with RHOEND<=RHOBEG. !! Typically, RHOBEG should be about one tenth of the greatest expected !! change to a variable, and RHOEND should indicate the accuracy that !! is required in the final values of the variables. real(wp),intent(in) :: rhoend !! RHOBEG and RHOEND must be set to the initial and final values of a !! trust region radius, so both must be positive with RHOEND<=RHOBEG. !! Typically, RHOBEG should be about one tenth of the greatest expected !! change to a variable, and RHOEND should indicate the accuracy that !! is required in the final values of the variables. integer,intent(in) :: iprint !! The value of IPRINT should be set to 0, 1, 2 or 3, which controls the !! amount of printing. Specifically, there is no output if IPRINT=0 and !! there is output only at the return if IPRINT=1. Otherwise, the best !! feasible vector of variables so far and the corresponding value of !! the objective function are printed whenever RHO is reduced, where !! RHO is the current lower bound on the trust region radius. Further, !! each new value of F with its variables are output if IPRINT=3. integer,intent(in) :: maxfun !! an upper bound on the number of calls of CALFUN, !! its value being at least NPT+1. procedure (func) :: calfun !! It must set !! F to the value of the objective function for the variables X(1), !! X(2),...,X(N). The value of the argument F is positive when CALFUN !! is called if and only if the current X satisfies the constraints real(wp),parameter :: zero = 0.0_wp real(wp),dimension(:),allocatable :: w integer, dimension(n) :: iact !to avoid type mismatch error - JW real(wp) :: smallx,sum,temp integer :: np,nptm,iamat,ib,iflag,i,iac,ibmat,ifv,igo,ihq,ipq,ipqw,& iqf,irc,isp,istp,iw,ixb,ixn,ixo,ixp,ixs,irf,izmat,j,ndim ! W is an array used for working space. Its length must be at least ! M*(2+N) + NPT*(4+N+NPT) + N*(9+3*N) + MAX [ M+3*N, 2*M+N, 2*NPT ]. ! On return, W(1) is set to the final value of F, and W(2) is set to ! the total number of function evaluations plus 0.5. allocate(w(M*(2+N) + NPT*(4+N+NPT) + N*(9+3*N) + MAX(M+3*N, 2*M+N, 2*NPT))) ! ! Check that N, NPT and MAXFUN are acceptable. ! smallx = 1.0e-6_wp * rhoend np = n + 1 nptm = npt - np if (n <= 1) then print 10 10 format (/ 4 x, 'Return from LINCOA because N is less than 2.') return end if if (npt < n+2 .or. npt > ((n+2)*np)/2) then print 20 20 format (/ 4 x, 'Return from LINCOA because NPT is not in',& ' the required interval.') return end if if (maxfun <= npt) then print 30 30 format (/ 4 x, 'Return from LINCOA because MAXFUN is less', ' than NPT+1.') return end if ! ! Normalize the constraints, and copy the resultant constraint matrix ! and right hand sides into working space, after increasing the right ! hand sides if necessary so that the starting point is feasible. ! iamat = max (m+3*n, 2*m+n, 2*npt) + 1 ib = iamat + m * n iflag = 0 if (m > 0) then iw = iamat - 1 do j = 1, m sum = zero temp = zero do i = 1, n sum = sum + a (i, j) * x (i) temp = temp + a (i, j) ** 2 end do if (temp == zero) then print 50 50 format (/ 4 x, 'Return from LINCOA because the gradient of',& ' a constraint is zero.') return end if temp = sqrt (temp) if (sum-b(j) > smallx*temp) iflag = 1 w (ib+j-1) = max (b(j), sum) / temp do i = 1, n iw = iw + 1 w (iw) = a (i, j) / temp end do end do end if if (iflag == 1) then if (iprint > 0) print 70 70 format (/ 4 x, 'LINCOA has made the initial X feasible by',& ' increasing part(s) of B.') end if ! ! Partition the working space array, so that different parts of it can be ! treated separately by the subroutine that performs the main calculation. ! ndim = npt + n ixb = ib + m ixp = ixb + n ifv = ixp + n * npt ixs = ifv + npt ixo = ixs + n igo = ixo + n ihq = igo + n ipq = ihq + (n*np) / 2 ibmat = ipq + npt izmat = ibmat + ndim * n istp = izmat + npt * nptm isp = istp + n ixn = isp + npt + npt iac = ixn + n irc = iac + n iqf = irc + m irf = iqf + n * n ipqw = irf + (n*np) / 2 ! ! The above settings provide a partition of W for subroutine LINCOB. ! call lincob (n, npt, m, w(iamat), w(ib), x, rhobeg, rhoend, iprint, maxfun, & w(ixb), w(ixp), w(ifv), w(ixs), w(ixo), w(igo), w(ihq), w(ipq), & w(ibmat), w(izmat), ndim, w(istp), w(isp), w(ixn), iact, w(irc), & w(iqf), w(irf), w(ipqw), w, calfun) deallocate(w) end subroutine lincoa