This subroutine minimizes an objective function F(X) subject to M inequality constraints on X, where X is a vector of variables that has N components. The algorithm employs linear approximations to the objective and constraint functions, the approximations being formed by linear interpolation at N+1 points in the space of the variables. We regard these interpolation points as vertices of a simplex. The parameter RHO controls the size of the simplex and it is reduced automatically from RHOBEG to RHOEND. For each RHO the subroutine tries to achieve a good vector of variables for the current size, and then RHO is reduced until the value RHOEND is reached. Therefore RHOBEG and RHOEND should be set to reasonable initial changes to and the required accuracy in the variables respectively, but this accuracy should be viewed as a subject for experimentation because it is not guaranteed.
The subroutine has an advantage over many of its competitors, however, which is that it treats each constraint individually when calculating a change to the variables, instead of lumping the constraints together into a single penalty function.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | n |
number of variables |
||
| integer, | intent(in) | :: | m |
number of inequality constraints |
||
| real(kind=wp), | intent(inout), | dimension(*) | :: | x |
Initial values of the variables must be set in X(1),X(2),...,X(N). On return they will be changed to the solution. |
|
| real(kind=wp), | intent(in) | :: | rhobeg |
reasonable initial change to variables (see description of RHO) |
||
| real(kind=wp), | intent(in) | :: | rhoend |
required accuracy (see description of RHO) |
||
| integer, | intent(in) | :: | iprint |
IPRINT should be set to 0, 1, 2 or 3, which controls the amount of printing during the calculation. Specifically, there is no output if IPRINT=0 and there is output only at the end of the calculation if IPRINT=1. Otherwise each new value of RHO and SIGMA is printed. Further, the vector of variables and some function information are given either when RHO is reduced or when each new value of F(X) is computed in the cases IPRINT=2 or IPRINT=3 respectively. Here SIGMA is a penalty parameter, it being assumed that a change to X is an improvement if it reduces the merit function F(X)+SIGMA*MAX(0.0,-C1(X),-C2(X),...,-CM(X)), where C1,C2,...,CM denote the constraint functions that should become nonnegative eventually, at least to the precision of RHOEND. In the printed output the displayed term that is multiplied by SIGMA is called MAXCV, which stands for 'MAXimum Constraint Violation'. |
||
| integer, | intent(inout) | :: | maxfun |
MAXFUN is an integer variable that must be set by the user to a limit on the number of calls of CALCFC. The value of MAXFUN will be altered to the number of calls of CALCFC that are made. |
||
| procedure(func) | :: | calcfc |
In order to define the objective and constraint functions, we require a subroutine that has the name and arguments SUBROUTINE CALCFC (N,M,X,F,CON) DIMENSION X(),CON() The values of N and M are fixed and have been defined already, while X is now the current vector of variables. The subroutine should return the objective and constraint functions at X in F and CON(1),CON(2), ...,CON(M). Note that we are trying to adjust X so that F(X) is as small as possible subject to the constraint functions being nonnegative. |
subroutine cobyla (n, m, x, rhobeg, rhoend, iprint, maxfun, calcfc) implicit none integer,intent(in) :: n !! number of variables integer,intent(in) :: m !! number of inequality constraints real(wp),dimension(*),intent(inout) :: x !! Initial values of the variables must be set in X(1),X(2),...,X(N). !! On return they will be changed to the solution. real(wp),intent(in) :: rhobeg !! reasonable initial change to variables (see description of RHO) real(wp),intent(in) :: rhoend !! required accuracy (see description of RHO) integer,intent(in) :: iprint !! IPRINT should be set to 0, 1, 2 or 3, which controls the amount of !! printing during the calculation. Specifically, there is no output if !! IPRINT=0 and there is output only at the end of the calculation if !! IPRINT=1. Otherwise each new value of RHO and SIGMA is printed. !! Further, the vector of variables and some function information are !! given either when RHO is reduced or when each new value of F(X) is !! computed in the cases IPRINT=2 or IPRINT=3 respectively. Here SIGMA !! is a penalty parameter, it being assumed that a change to X is an !! improvement if it reduces the merit function !! F(X)+SIGMA*MAX(0.0,-C1(X),-C2(X),...,-CM(X)), !! where C1,C2,...,CM denote the constraint functions that should become !! nonnegative eventually, at least to the precision of RHOEND. In the !! printed output the displayed term that is multiplied by SIGMA is !! called MAXCV, which stands for 'MAXimum Constraint Violation'. integer,intent(inout) :: maxfun !! MAXFUN is an integer variable that must be set by the user to a !! limit on the number of calls of CALCFC. !! The value of MAXFUN will be altered to the number of calls !! of CALCFC that are made. procedure (func) :: calcfc !! In order to define the objective and constraint functions, we require !! a subroutine that has the name and arguments !! SUBROUTINE CALCFC (N,M,X,F,CON) !! DIMENSION X(*),CON(*) !! The values of N and M are fixed and have been defined already, while !! X is now the current vector of variables. The subroutine should return !! the objective and constraint functions at X in F and CON(1),CON(2), !! ...,CON(M). Note that we are trying to adjust X so that F(X) is as !! small as possible subject to the constraint functions being nonnegative. integer,dimension(:),allocatable :: iact real(wp),dimension(:),allocatable :: w integer :: mpp,icon,isim,isimi,idatm,ia,ivsig,iveta,isigb,idx,iwork !W and IACT provide real and integer arrays that are used as working space. allocate(w(N*(3*N+2*M+11)+4*M+6)) allocate(iact(M+1)) ! Partition the working space array W to provide the storage that is needed ! for the main calculation. mpp = m + 2 icon = 1 isim = icon + mpp isimi = isim + n * n + n idatm = isimi + n * n ia = idatm + n * mpp + mpp ivsig = ia + m * n + n iveta = ivsig + n isigb = iveta + n idx = isigb + n iwork = idx + n call cobylb (n, m, mpp, x, rhobeg, rhoend, iprint, maxfun, w(icon), w(isim), & w(isimi), w(idatm), w(ia), w(ivsig), w(iveta), w(isigb), w(idx), & w(iwork), iact, calcfc) deallocate(iact) deallocate(w) end subroutine cobyla