True BASIC Program to Find a Function Minimum

Using the Powell Search Method

by Namir Shammas

The following program uses the Powell Search method to find the minimum of a function.

The program prompts you to enter for each variable (i.e. dimension):

1. The tolerance for the minimized function,

2. The tolerance for the search step.

3. The maximum number of iterations.

4. The initial guesses for the optimum point for each variable.

The program displays intermediate values for the function and the variables. The program displays the following final results:

1. The coordinates of the minimum point.

2. The minimum function value.

3. The number of iterations

The current code finds the minimum for the following function:

f(x1,x2) = x1 - x2 + 2 * x1 ^ 2 + 2 * x1 * x2 + x2 ^ 2

Here is a sample session to solve for the optimum of the above function:

Here is the BASIC listing:

! Powell Search Optimization
OPTION TYPO
OPTION NOLET

DECLARE NUMERIC MAX_VARS
DECLARE NUMERIC N, I, J, bTrue, bFalse 
DECLARE NUMERIC K, L, bOK 
DECLARE NUMERIC MaxCycles 
DECLARE NUMERIC F1, F2, F 
DECLARE NUMERIC Eps_Fx, Eps_Step 
DECLARE NUMERIC Alpha, Sum 

DIM X(1), X1(1), X2(1) 
DIM S(1, 1)

MAX_VARS = 2
bTrue = 1
bFalse = 0

SUB MyFx(X(), N, Res)
  ! Res = 100 * (X(1) ^ 2 - X(2)) ^ 2 + (1 - X(1)) ^ 2
  Res = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2
End SUB

SUB MyFxEx(N, X(), S(,), II, Alpha, funRes)
  LOCAL I
  DIM XX(1) 
  
  MAT REDIM XX(N)
  
  For I = 1 To N
    XX(I) = X(I) + Alpha * S(II, I)
  Next I
  
  CALL MyFx(XX, N, funRes)
End SUB

Sub GetGradients(X(), N, Deriv(), DerivNorm)
  LOCAL XX, I, H, Fp, Fm

  DerivNorm = 0
  For I = 1 To N
    XX = X(I)
    H = 0.01 * (1 + ABS(XX))
    X(I) = XX + H
    CALL MyFx(X, N, Fp)
    X(I) = XX - H
    CALL MyFx(X, N, Fm)
    X(I) = XX
    Deriv(I) = (Fp - Fm) / 2 / H
    DerivNorm = DerivNorm + Deriv(I) ^ 2
  Next I
  DerivNorm = Sqr(DerivNorm)
End Sub

SUB LinSearch_DirectSearch(X(), N, Alpha, S(,), II, InitStep, MinStep, boolRes)
  LOCAL F1, F2
  
  CALL MyFxEx(N, X, S, II, Alpha, F1)
  
  Do
    CALL MyFxEx(N, X, S, II, Alpha + InitStep, F2)
    If F2 < F1 Then
      F1 = F2
      Alpha = Alpha + InitStep
    Else
      CALL MyFxEx(N, X, S, II, Alpha - InitStep, F2)
      If F2 < F1 Then
        F1 = F2
        Alpha = Alpha - InitStep
      Else
        ! reduce search step size
        InitStep = InitStep / 10
      End If
    End If
  Loop Until InitStep < MinStep
  
  boolRes = bTrue

End SUB

! Powell Search Optimization

MAT REDIM X(MAX_VARS) 
MAT REDIM X1(MAX_VARS) 
MAT REDIM X2(MAX_VARS) 
MAT REDIM S(MAX_VARS+1, MAX_VARS) 

N = MAX_VARS

PRINT "Powell Search Optimization"
INPUT PROMPT "Enter function tolerance  value: ": Eps_Fx
INPUT PROMPT "Enter step tolerance  value: ": Eps_Step
INPUT PROMPT "Enter maximum number of cycles: ": MaxCycles
For I = 1 To N
  PRINT "Enter guess for X(";I;")";
  INPUT X(I)
Next I

J = 1
CALL MyFx(X, N, F1)
For K = 1 To N
  X1(K) = X(K)
  For L = 1 To N
    S(K,L) = 0
    IF K=L Then S(K, L) = 1
  Next L
Next K

Do
  
  ! reset row S(N+1,:)
  For K = 1 To N
    S(N + 1, K) = 0
  Next K
  
  For I = 1 To N
    CALL LinSearch_DirectSearch(X, N, Alpha, S, I, 0.1, 0.0001, bOK)
    For K = 1 To N
      X(K) = X(K) + Alpha * S(I, K)
      S(N + 1, K) = S(N + 1, K) + Alpha * S(I, K)
    Next K
  Next I
  
  CALL LinSearch_DirectSearch(X, N, Alpha, S, N + 1, 0.1, 0.0001, bOK)
  
  For K = 1 To N
    X(K) = X(K) + Alpha * S(N + 1, K)
    X2(K) = X(K)
  Next K
  CALL MyFx(X2, N, F2)

  PRINT "F = ";F2;" ";
  For K = 1 To N
    PRINT "X(";K;")=";X(K);" ";
  Next K
  PRINT
  
  ! test end of iterations criteria
  If Abs(F2 - F1) < Eps_Fx Then Exit Do
  Sum = 0
  For K = 1 To N
    Sum = Sum + (X2(K) - X1(K)) ^ 2
  Next K
  Sum = Sqr(Sum)
  If Sum < Eps_Step Then Exit Do
  If J => MaxCycles Then Exit Do
  
  J = J + 1
  
  ! rotate data
  For K = 1 To N
    X1(K) = X2(K)
  Next K
  F1 = F2
  
  ! copy S matrix
  For K = 1 To N
    For L = 1 To N
      S(K, L) = S(K + 1, L)
    Next L
  Next K

Loop

PRINT "**********FINAL RESULTS************"
PRINT "Optimum at:"
For I = 1 To N
  PRINT "X(";I;")=";X(I)
Next I
PRINT "Function value ="; F1
PRINT "Number of iterations = ";J

END

BACK