cobyla.f 4.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475
  1. SUBROUTINE COBYLA (N,M,X,RHOBEG,RHOEND,IPRINT,MAXFUN,W,IACT,ierr)
  2. IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  3. DIMENSION X(*),W(*),IACT(*)
  4. C
  5. C This subroutine minimizes an objective function F(X) subject to M
  6. C inequality constraints on X, where X is a vector of variables that has
  7. C N components. The algorithm employs linear approximations to the
  8. C objective and constraint functions, the approximations being formed by
  9. C linear interpolation at N+1 points in the space of the variables.
  10. C We regard these interpolation points as vertices of a simplex. The
  11. C parameter RHO controls the size of the simplex and it is reduced
  12. C automatically from RHOBEG to RHOEND. For each RHO the subroutine tries
  13. C to achieve a good vector of variables for the current size, and then
  14. C RHO is reduced until the value RHOEND is reached. Therefore RHOBEG and
  15. C RHOEND should be set to reasonable initial changes to and the required
  16. C accuracy in the variables respectively, but this accuracy should be
  17. C viewed as a subject for experimentation because it is not guaranteed.
  18. C The subroutine has an advantage over many of its competitors, however,
  19. C which is that it treats each constraint individually when calculating
  20. C a change to the variables, instead of lumping the constraints together
  21. C into a single penalty function. The name of the subroutine is derived
  22. C from the phrase Constrained Optimization BY Linear Approximations.
  23. C
  24. C The user must set the values of N, M, RHOBEG and RHOEND, and must
  25. C provide an initial vector of variables in X. Further, the value of
  26. C IPRINT should be set to 0, 1, 2 or 3, which controls the amount of
  27. C printing during the calculation. Specifically, there is no output if
  28. C IPRINT=0 and there is output only at the end of the calculation if
  29. C IPRINT=1. Otherwise each new value of RHO and SIGMA is printed.
  30. C Further, the vector of variables and some function information are
  31. C given either when RHO is reduced or when each new value of F(X) is
  32. C computed in the cases IPRINT=2 or IPRINT=3 respectively. Here SIGMA
  33. C is a penalty parameter, it being assumed that a change to X is an
  34. C improvement if it reduces the merit function
  35. C F(X)+SIGMA*MAX(0.0,-C1(X),-C2(X),...,-CM(X)),
  36. C where C1,C2,...,CM denote the constraint functions that should become
  37. C nonnegative eventually, at least to the precision of RHOEND. In the
  38. C printed output the displayed term that is multiplied by SIGMA is
  39. C called MAXCV, which stands for 'MAXimum Constraint Violation'. The
  40. C argument MAXFUN is an integer variable that must be set by the user to a
  41. C limit on the number of calls of CALCFC, the purpose of this routine being
  42. C given below. The value of MAXFUN will be altered to the number of calls
  43. C of CALCFC that are made. The arguments W and IACT provide real and
  44. C integer arrays that are used as working space. Their lengths must be at
  45. C least N*(3*N+2*M+11)+4*M+6 and M+1 respectively.
  46. C
  47. C In order to define the objective and constraint functions, we require
  48. C a subroutine that has the name and arguments
  49. C SUBROUTINE CALCFC (N,M,X,F,CON)
  50. C DIMENSION X(*),CON(*) .
  51. C The values of N and M are fixed and have been defined already, while
  52. C X is now the current vector of variables. The subroutine should return
  53. C the objective and constraint functions at X in F and CON(1),CON(2),
  54. C ...,CON(M). Note that we are trying to adjust X so that F(X) is as
  55. C small as possible subject to the constraint functions being nonnegative.
  56. C
  57. C Partition the working space array W to provide the storage that is needed
  58. C for the main calculation.
  59. C
  60. MPP=M+2
  61. ICON=1
  62. ISIM=ICON+MPP
  63. ISIMI=ISIM+N*N+N
  64. IDATM=ISIMI+N*N
  65. IA=IDATM+N*MPP+MPP
  66. IVSIG=IA+M*N+N
  67. IVETA=IVSIG+N
  68. ISIGB=IVETA+N
  69. IDX=ISIGB+N
  70. IWORK=IDX+N
  71. CALL COBYLB (N,M,MPP,X,RHOBEG,RHOEND,IPRINT,MAXFUN,W(ICON),
  72. 1 W(ISIM),W(ISIMI),W(IDATM),W(IA),W(IVSIG),W(IVETA),W(ISIGB),
  73. 2 W(IDX),W(IWORK),IACT, ierr)
  74. RETURN
  75. END