minimisation - Miminisation routines

The *newuoa_minimisation* and *bobyqa_minimisation* modules provide facilities
to minimise a function of one or more variables without the use of derivatives
(see A view of algorithms for optimization without derivatives.

The first module provides the routine *NEWUOA* which implements an unconstrained minimisation
algorithm. The second module provides the routine *BOBYQA* which implements an algorithm where
the minimum of the function is sought within a given hyperblock.

*Note:* The code was developed by Mike Powell. It has been slightly adjusted so that advantage
can be taken of a number of modern Fortran features. This has led to a slightly different interface
than the original code. These are the changes:

Use the generic names for standard functions like

*abs*and*min*.Use an allocatable work array, so that the user is no longer responsible for this.

Wrap the code in a module, so that the interfaces can be checked by the compiler.

Use a

*kind*parameter instead of*REAL*8*.Use assumed shape arrays, so that the number of dimensions can be automatically derived from the size of the

*X*array.Make the subroutine that computes the value of the function to be minimised an argument, rather than requiring a fixed name. This makes the module more flexible.

There are two public routines, one in each module:

**call newuoa((npt,x,rhobeg,rhoend,iprint,maxfun,calfun)**This subroutine seeks the least value of a function of many variables, by a trust region method that forms quadratic models by interpolation. There can be some freedom in the interpolation conditions, which is taken up by minimizing the Frobenius norm of the change to the second derivative of the quadratic model, beginning with a zero matrix. The arguments of the subroutine are as follows.

Initial values of the variables must be set in

*X(1),X(2),...,X(N)*. They will be changed to the values that give the least calculated*F*.- integer
*npt* *NPT*is the number of interpolation conditions. Its value must be in the interval [*N+2,(N+1)(N+2)/2*].- real(kind=wp), dimension(:)
*x* Array with coordinate values. On input the initial values, on output the vector for which the minimum function value was found.

- real(kind=wp)
*rhobeg, 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
*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, each new value of*RHO*is printed, with the best vector of variables so far and the corresponding value of the objective function. Further, each new value of*F*with its variables are output if*IPRINT=3*.- integer
*maxfun* *MAXFUN*must be set to an upper bound on the number of calls of*CALFUN*.- subroutine
*calfun* Subroutine

*CALFUN*must be provided by the user. It must set*F*to the value of the objective function for the variables*X(1),X(2),...,X(N)*.Its interface is:

subroutine calfun( x, f ) import wp real(kind=wp), dimension(:) :: x real(kind=wp) :: f end subroutine calfun

- integer
**call bobyqa((npt,x,xl,xu,rhobeg,rhoend,iprint,maxfun,calfun)**This subroutine seeks the least value of a function of many variables, by applying a trust region method that forms quadratic models by interpolation. There is usually some freedom in the interpolation conditions, which is taken up by minimizing the Frobenius norm of the change to the second derivative of the model, beginning with the zero matrix. The values of the variables are constrained by upper and lower bounds. The arguments of the subroutine are as follows.

Initial values of the variables must be set in

*X(1),X(2),...,X(N)*. They will be changed to the values that give the least calculated*F*.- integer
*npt* *NPT*is the number of interpolation conditions. Its value must be in the interval [*N+2,(N+1)(N+2)/2*].- real(kind=wp), dimension(:)
*x* Array with coordinate values. On input the initial values, on output the vector for which the minimum function value was found.

- real(kind=wp), dimension(:)
*xl, xu* For

*I=1,2,...,N*,*XL(I)*and*XU(I)*must provide the lower and upper bounds, respectively, on*X(I)*. The construction of quadratic models requires*XL(I)*to be strictly less than*XU(I)*for each*I*. Further, the contribution to a model from changes to the*I*-th variable is damaged severely by rounding errors if*XU(I)-XL(I)*is too small.- real(kind=wp)
*rhobeg, 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*no greater than*RHOBEG*. Typically,*RHOBEG*should be about one tenth of the greatest expected change to a variable, while*RHOEND*should indicate the accuracy that is required in the final values of the variables. An error return occurs if any of the differences*XU(I)-XL(I), I=1,...,N*is less than*2*RHOBEG*.- integer
*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, each new value of*RHO*is printed, with the best vector of variables so far and the corresponding value of the objective function. Further, each new value of*F*with its variables are output if*IPRINT=3*.- integer
*maxfun* *MAXFUN*must be set to an upper bound on the number of calls of*CALFUN*.- subroutine
*calfun* Subroutine

*CALFUN*must be provided by the user. It must set*F*to the value of the objective function for the variables*X(1),X(2),...,X(N)*.Its interface is the same for routine

*NEWUOA*.

- integer

The functon to be minimised does not need to be a straightforward mathematical function. The only requirements are that it is continuous in its arguments and more or less smooth.

Here is an example of a minimisation problem that would be very difficult to solve when
derivatives are required. The problem is to find the values for the coefficients *a* and
*b* in the system of ordinary differential equations below that lead to a final vector
*(x,y) = (0.5,0.25)* at time *t = 1*:

dx/dt = - a x dy/dt = + a x - b y (x,y) = (1,0) at time t = 0

The subroutine to compute the function to be minimised is shown below. It computes the
solution at *t = 1* for the above system and computes the sum of square differences with the
target vector. This sum is the function value to be minimised:

! compute_diff -- ! Compute the solution and determine how close it is to the required solution ! ! Arguments: ! param Vector of parameters ! f Value of the function to be minimized ! subroutine compute_diff( param, f ) use bobyqa_minimisation, only: wp implicit none real(kind=wp), dimension(:) :: param real(kind=wp) :: f integer :: i real(kind=wp), dimension(2) :: x real(kind=wp), dimension(2) :: dx real(kind=wp) :: dt ! ! Initial condition: x = [1.0, 0.0] ! We want x = [0.5, 0.25] at time = 1 ! dt = 0.01 x(1) = 1.0 x(2) = 0.0 do i = 1,100 dx(1) = -param(1) * x(1) dx(2) = +param(1) * x(1) - param(2) * x(2) x(1) = x(1) + dx(1) * dt x(2) = x(2) + dx(2) * dt enddo f = (x(1)-0.5)**2 + (x(2)-0.25)**2 end subroutine calfun

Copyright © 2013 Mike Powell

Copyright © 2013 Arjen Markus <arjenmarkus at sourceforge dot net>