/*************************************************************************/ /* */ /* ALGORITHM AS 298: A HYBRID OPTIMIZATION ALGORITHM */ /* */ /* Stephen P. Brooks */ /* University of Kent, Canterbury, England. */ /* */ /*************************************************************************/ /* Journal of the Royal Statistical Society (Series C): */ /* Applied Statistics, Vol.44, No. 4, 1995, pp. 530-533. */ /* Translated from Fortran to PL/I by R. A. Vowels, April 2006. */ (SUBSCRIPTRANGE, FIXEDOVERFLOW, SIZE): hybrid: PROCEDURE (nparam, param, n, t0, rho, nt, bl, bu, np, points, low, ifault) OPTIONS (REORDER); /* .. Scalar Arguments .. */ DECLARE ( nparam ) FIXED BINARY (31); DECLARE ( param(*) ) FLOAT (18); DECLARE ( n ) FIXED BINARY (31); DECLARE ( t0 ) FLOAT (18); DECLARE ( rho ) FLOAT (18); DECLARE ( nt ) FIXED BINARY (31); DECLARE ( bl(*) ) FLOAT (18); DECLARE ( bu(*) ) FLOAT (18); DECLARE ( np ) FIXED BINARY (31); DECLARE ( points(*,*) ) FLOAT (18); /* points(nparam,n) */ DECLARE ( low(*) ) FLOAT (18); DECLARE ( ifault ) FIXED BINARY (31); /* .. */ /* .. Local Scalars .. */ DECLARE ( bst, diff, newobj, objf, p, t, u ) FLOAT (18); DECLARE ( i,i1,i2 ) FIXED BINARY (31); /* .. */ /* Define the constants ZERO and ONE and set */ /* initial values for BST and IFAULT */ /* .. Data statements .. */ DECLARE ( zero STATIC INITIAL ( 0.00000000000000000E0), one STATIC INITIAL ( 1.00000000000000000E0) ) FLOAT (18); /* .. */ ifault = 0; bst = 1.00000000000000000E+10; /* ****************************************** */ /* *** Start of Annealing component to *** */ /* *** generate a set of starting points. *** */ /* ****************************************** */ /* Initialise the temperature */ t = t0; /* Choose an initial point */ DO i = 1 TO nparam; param(i) = random(bl(i), bu(i)); END; /* Calculate the function value for the initial point */ CALL funct1(nparam, param, objf); /* Begin the loop over I1 until Nt Temp. reductions. */ DO i1 = 1 TO nt; /* Set the counter for the No of points accepted to one */ np = 1; /* Begin the loop over I2 until N points have been considered. */ DO i2 = 1 TO n; /* Record the present point */ DO i = 1 TO nparam; points(i,np+1) = param(i); END; /* Select a new point */ DO i = 1 TO nparam; param(i) = random(bl(i), bu(i)); END; /* Calculate the function value of this new point */ CALL funct1(nparam, param, newobj); /* Calculate the difference */ diff = newobj - objf; /* If the new point is better than the old one then increase the */ /* Np counter, and record it as a possible starting point. Move */ /* to this new point. */ IF diff < zero THEN DO; np = np + 1; DO i = 1 TO nparam; points(i,np) = param(i); END; objf = newobj; IF newobj < bst THEN DO; bst = newobj; DO i = 1 TO nparam; points(i,1) = param(i); END; END; /* Otherwise use the Metropolis Criterion */ END; ELSE DO; /* P is the comparison statistic. U is a standard Uniform variate. */ p = EXP(-diff/t); u = random(zero, one); /* If U < P then accept this point. Increase the Np counter. */ /* Record this point as a possible starting point, and move */ /* to this new point. */ IF u < p THEN DO; np = np + 1; objf = newobj; DO i = 1 TO nparam; points(i,np) = param(i); END; /* If U >= P then reject this point and return to the old point. */ END; ELSE DO; DO i = 1 TO nparam; param(i) = points(i,np+1); END; END; /* Both alternatives for the comparison of the new and old point */ /* Have been considered. */ END; /* Continue with the loop until N new points have been considered */ END; /* Set the error indicator to 1 if the algorithm appears to have converged. */ IF np = 1 THEN ifault = i1; /* Once N points have been considered, lower the Temperature. */ t = t*rho; /* Continue with the loop until the Temp. has been reduced Nt times */ END; /* ***************************************** */ /* *** End of Annealing component, NP *** */ /* *** starting points are now stored in *** */ /* *** the array POINTS. *** */ /* ***************************************** */ /* Set an initially high value of LOW(NPARAM+1) */ low(nparam+1) = 1.00000000000000000E+10; /* ******************************************* */ /* *** Start of the traditional component, *** */ /* *** which loops over the NP starting *** */ /* *** points calling the traditional *** */ /* *** minimisation routine for each. *** */ /* ******************************************* */ /* Begin the loop over the Np starting points. */ DO i1 = 1 TO np; /* Set the PARAM array to be the values of the I1st start point */ DO i2 = 1 TO nparam; param(i2) = points(i2,i1); END; /* Call traditional minimisation routine */ CALL tradnl(nparam, param, bl, bu, objf); /* Check to see if this starting point gives the */ /* lowest function value to date */ IF objf < low(nparam+1) THEN DO; DO i = 1 TO nparam; low(i) = param(i); END; low(nparam+1) = objf; END; /* End the Np loop */ END; /* ************************************** */ /* *** End of the second, traditional *** */ /* *** component. *** */ /* ************************************** */ /* Return to calling routine */ RETURN; /* End of subroutine HYBRID */ random: PROCEDURE (bl, bu) RETURNS( FLOAT (18) ) OPTIONS (REORDER); DECLARE ( bl, bu ) FLOAT (18); DECLARE RANDOM BUILTIN; RETURN ( bl + (bu - bl)*RANDOM ); END random; END hybrid;