'************************************************************** '* Minimization of a Function FUNC of N Variables By Powell's * '* Method Discarding the Direction of Largest Decrease * '* ---------------------------------------------------------- * '* SAMPLE RUN: Find a minimum of function F(x,y): * '* F=Sin(R)/R, where R = Sqr(x*x+y*y). * '* * '* Number of iterations: 2 * '* * '* Minimum Value: -.2172336282112217 * '* * '* at point: 3.177320295222616 3.177320295222616 * '* * '* ---------------------------------------------------------- * '* REFERENCE: "Numerical Recipes, The Art of Scientific * '* Computing By W.H. Press, B.P. Flannery, * '* S.A. Teukolsky and W.T. Vetterling, * '* Cambridge University Press, 1986" * '* [BIBLI 08]. * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '************************************************************** 'PROGRAM TPOWELL DEFDBL A-H, O-Z DEFINT I-N NP = 20 DIM P(NP), XI(NP, NP), P1(NP) DIM PCOM(NP), XICOM(NP) 'PCOM,XICOM,NCOM are common variables DIM XT(NP) 'for LINMIN and F1DIM only N = 2 P(1) = 2#: P(2) = 2# XI(1, 1) = 1#: XI(1, 2) = 1# XI(2, 1) = 1#: XI(2, 2) = 1# FTOL = .000001# GOSUB 1000 'call POWELL(P,XI,N,NP,FTOL,ITER,FRET) CLS PRINT PRINT " Number of iterations: "; ITER PRINT PRINT " Minimum value: "; FRET PRINT PRINT " at point: "; P(1); " "; P(2) PRINT END 'of main program 'user defined function to minimize 500 'FUNCTION FUNC(P1) R1 = SQR(P1(1) * P1(1) + P1(2) * P1(2)) IF ABS(R1) < .000000000001# THEN FUNC = 1# ELSE FUNC = SIN(R1) / R1 END IF RETURN 600 'FUNCTION XMAX(a,b) IF A >= B THEN XMAX = A ELSE XMAX = B END IF RETURN 610 'FUNCTION XMIN(a,b) IF A <= B THEN XMIN = A ELSE XMIN = B END IF RETURN 650 'FUNCTION Sign(a,b) IF B >= 0 THEN Sign = ABS(A) ELSE Sign = -ABS(A) END IF RETURN 1000 'Subroutine POWELL(P, XI, N, FTOL, ITER, FRET) '----------------------------------------------------------- ' Minimization of a function FUNC of N variables (FUNC is ' not an argument, it is a fixed function name). Input con- ' sists of an initial starting point P that is a vector of ' length N: an initial matrix XI whose logical dimensions ' are N by N, physical dimensions NP by NP, and whose columns ' contain the initial set of directions (usually the N unit ' vectors): and FTOL, the fractional tolerance in the func- ' tion value such that failure to decrease by more than this ' amount on one iteration signals doneness. On output, P is ' set to the best point found, XI is the then-current direc- ' tion set, FRET is the returned function value at P, and ' ITER is the number of iterations taken. The routine LINMIN ' is used. '----------------------------------------------------------- ' Label: 1 ITMAX = 200 DIM PT(NP), PTT(NP), XIT(NP) 'FRET=FUNC(P) P1(1) = P(1): P1(2) = P(2): GOSUB 500: FRET = FUNC FOR J = 1 TO N PT(J) = P(J): 'Save initial point NEXT J ITER = 0 1 ITER = ITER + 1 FP = FRET IBIG = 0 DEL = 0# 'Will be the biggest function decrease. FOR I = 1 TO N 'In each iteration, loop over all directions in the set. FOR J = 1 TO N 'Copy the direction. XIT(J) = XI(J, I) NEXT J FPTT = FRET GOSUB 1500 'call LINMIN(P,XIT,N,FRET) (Minimize along it) IF ABS(FPTT - FRET) > DEL THEN DEL = ABS(FPTT - FRET) IBIG = I END IF NEXT I IF 2# * ABS(FP - FRET) <= FTOL * (ABS(FP) + ABS(FRET)) THEN RETURN'Termination criterion IF ITER = ITMAX THEN PRINT " Powell exceeding maximum iterations." RETURN END IF FOR J = 1 TO N PTT(J) = 2# * P(J) - PT(J)'Construct the extrapolated point and the average XIT(J) = P(J) - PT(J) 'direction moved. Save the old starting point PT(J) = P(J) NEXT J 'FPTT=FUNC(PTT) (Function value at extrapolated point) P1(1) = PTT(1): P1(2) = PTT(2): GOSUB 500: FPTT = FUNC IF FPTT >= FP THEN GOTO 1: 'One reason not to use new direction. T = 2# * (FP - 2# * FRET + FPTT) * (FP - FRET - DEL) * (FP - FRET - DEL) - DEL * (FP - FPTT) * (FP - FPTT) IF T >= 0! THEN GOTO 1: 'Other reason not to use new direction. GOSUB 1500 'call LINMIN(P,XIT,N,FRET) 'Move to the minimum of the new direction. FOR J = 1 TO N 'and save the new direction. XI(J, IBIG) = XIT(J) NEXT J GOTO 1 RETURN 1500 'Subroutine LINMIN(P, XIT, N, FRET) '---------------------------------------------------------- ' Given an N dimensional point P and a N dimensional direc- ' tion XI, moves and resets P to where the function FUNC(P) ' takes on a minimum along the direction XIT from P, and ' replaces XIT by the actual vector displacement that P was ' moved. Also returns as FRET the value of FUNC at the ' returned location P. This is actually all accomplished by ' calling the routines MNBRAK and BRENT. '---------------------------------------------------------- TOL = .0001# NCOM = N FOR J = 1 TO N PCOM(J) = P(J) XICOM(J) = XIT(J) NEXT J AX = 0# BX = 1# CX = 2# GOSUB 2000 'call MNBRAK(AX,BX,CX,FA,FB,FC) GOSUB 3000 'call BRENT(AX,BX,CX,TOL,XMIN) FRET = BRENT FOR J = 1 TO N XIT(J) = XMIN * XIT(J) P(J) = P(J) + XIT(J) NEXT J RETURN 1600 'FUNCTION F1DIM(X1) FOR J = 1 TO NCOM XT(J) = PCOM(J) + X1 * XICOM(J) NEXT J 'F1DIM = FUNC(XT) P1(1) = XT(1): P1(2) = XT(2): GOSUB 500 F1DIM = FUNC RETURN 2000 'Subroutine MNBRAK(AX,BX,CX,FA,FB,FC) '-------------------------------------------------------------------- 'Given a Function F1DIM(X), and given distinct initial points AX and 'BX, this routine searches in the downhill direction (defined by the 'Function as evaluated at the initial points) and returns new points 'AX, BX, CX which bracket a minimum of the Function. Also returned 'are the Function values at the three points, FA, FB and FC. '-------------------------------------------------------------------- ' Label: 10 GOLD = 1.618034: GLIMIT = 100#: TINY = 9.999999999999999D-21 'The first parameter is the default ratio by which successive intervals 'are magnified: the second is the maximum magnification allowed for 'a parabolic-fit step. 'FA=F1DIM(AX) X1 = AX: GOSUB 1600: FA = F1DIM 'FB=F1DIM(BX) X1 = BX: GOSUB 1600: FB = F1DIM IF FB > FA THEN DUM = AX AX = BX BX = DUM DUM = FB FB = FA FA = DUM END IF CX = BX + GOLD * (BX - AX) 'FC=F1DIM(CX) X1 = CX: GOSUB 1600: FC = F1DIM 10 IF FB >= FC THEN R = (BX - AX) * (FB - FC) Q = (BX - CX) * (FB - FA) A = Q - R: B = TINY: GOSUB 600 A = XMAX: B = Q - R: GOSUB 650 U = BX - ((BX - CX) * Q - (BX - AX) * R) / (2# * Sign) ULIM = BX + GLIMIT * (CX - BX) IF (BX - U) * (U - CX) > 0! THEN 'FU=F1DIM(U) X1 = U: GOSUB 1600: FU = F1DIM IF FU < FC THEN AX = BX FA = FB BX = U FB = FU GOTO 10 ELSEIF FU > FB THEN CX = U FC = FU GOTO 10 END IF U = CX + GOLD * (CX - BX) 'FU=F1DIM(U) X1 = U: GOSUB 1600: FU = F1DIM ELSEIF (CX - U) * (U - ULIM) > 0! THEN 'FU=F1DIM(U) X1 = U: GOSUB 1600: FU = F1DIM IF FU < FC THEN BX = CX CX = U U = CX + GOLD * (CX - BX) FB = FC FC = FU 'FU=F1DIM(U) X1 = U: GOSUB 1600: FU = F1DIM END IF ELSEIF (U - ULIM) * (ULIM - CX) >= 0! THEN U = ULIM 'FU=F1DIM(U) X1 = U: GOSUB 1600: FU = F1DIM ELSE U = CX + GOLD * (CX - BX) 'FU=F1DIM(U) X1 = U: GOSUB 1600: FU = F1DIM END IF AX = BX BX = CX CX = U FA = FB FB = FC FC = FU GOTO 10 END IF RETURN 3000 'Function BRENTAX,BX,CX,TOL,XMIN) '------------------------------------------------------------------- 'Given a function F1DIM, and a bracketing triplet of abscissas 'AX,BX,CX (such that BX is between AX and CX, and F(BX) is less 'than both F(AX) and F(CX)), this routine isolates the minimum 'to a fractional precision of about TOL using Brent's method. 'The abscissa of the minimum is returned in XMIN, and the minimum 'function value is returned as BRENT, the returned function value. '------------------------------------------------------------------- 'Labels: 11,2,3 ITMAX = 100: CGOLD = .381966#: ZEPS = 1E-10 'Maximum allowed number of iterations: golden ratio: and a small 'number which protects against trying to achieve fractional accuracy 'for a minimum that happens to be exactly zero. 'A=MIN(AX,CX) IF AX <= CX THEN A = AX ELSE A = CX END IF 'B=MAX(AX,CX) IF AX >= CX THEN B = AX ELSE B = CX END IF V = BX W = V X = V E = 0# 'FX=F1DIM(X) X1 = X: GOSUB 1600: FX = F1DIM FV = FX FW = FX FOR ITER1 = 1 TO ITMAX 'main loop XM = .5 * (A + B) TOL1 = TOL * ABS(X) + ZEPS TOL2 = 2# * TOL1 IF ABS(X - XM) <= (TOL2 - .5 * (B - A)) THEN GOTO 3'Test for done here IF ABS(E) > TOL1 THEN 'Construct a trial parabolic fit R = (X - W) * (FX - FV) Q = (X - V) * (FX - FW) P = (X - V) * Q - (X - W) * R Q = .2 * (Q - R) IF Q > 0# THEN P = -P Q = ABS(Q) ETEMP = E E = D IF (ABS(P) >= ABS(.5 * Q * ETEMP)) OR (P <= Q * (A - X)) OR (P >= Q * (B - X)) THEN GOTO 11 ' The above conditions determine the acceptability of the ' parabolic fit. Here it is o.k.: D = P / Q U = X + D IF XM - X >= 0 THEN Sign = ABS(TOL1) ELSE Sign = -ABS(TOL1) END IF IF (U - A < TOL2) OR (B - U < TOL2) THEN D = Sign GOTO 2 END IF 11 IF X >= XM THEN E = A - X ELSE E = B - X END IF D = CGOLD * E 2 IF ABS(D) >= TOL1 THEN U = X + D ELSE IF D >= 0 THEN Sign = ABS(TOL1) ELSE Sign = -ABS(TOL1) END IF U = X + Sign END IF 'FU=F1DIM(U) (This the one function evaluation per iteration) X1 = U: GOSUB 1600: FU = F1DIM IF FU <= FX THEN IF U >= X THEN A = X ELSE B = X END IF V = W FV = FW W = X FW = FX X = U FX = FU ELSE IF U < X THEN A = U ELSE B = U END IF IF (FU <= FW) OR (W = X) THEN V = W FV = FW W = U FW = FU ELSEIF (FU <= FV) OR (V = X) OR (V = W) THEN V = U FV = FU END IF END IF NEXT ITER1 PRINT " Brent exceed maximum iterations." 3 XMIN = X 'exit section BRENT = FX RETURN 'end of file tpowell.bas