HP-71B Program to Find a Function Minimum

Using the Hooke-Jeeves Directional Search Method

by Namir Shammas

The following program calculates the minimum point of a multi-variable function using the Hooke-Jeeves directional search method.

The program prompts you to enter for each variable:

1. Guess for the minimum point.

2. Initial search step value.

3. The minimum search step value.

The program displays intermediate results to show the iteration progress and also the following final results:

1. The coordinates of the minimum value.

2. The step size for each variable.

3. The minimum function value.

4. 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. The first figure shows the input.:

The second figure shows the output:

Here is the BASIC listing:  

10 ! Directional Search using Hooke and Jeeves method
20 !
30 ! JPCROM-BASED VERSION 1.0
40 ! USES JPC ROM TO SUPPORT IF-ELSE-END IF AND LOOP-END LOOP STATEMENTS
50 !
60 !---------------------------- MAIN PROGRAM --------------------
70 !
80 N = 2 ! number of variables
90 DIM X(N), X1(N), M(N)
100 DIM S9(N), S0(N), D(N)
110 !
120 DISP "Hooke-Jeeves Directional Search Method"
130 For I = 1 To N
140 DISP "Enter value for X(";I;")";
150 INPUT X1(I)
160 DISP "Enter step size for X(";I;")";
170 INPUT S9(I)
180 DISP "Enter minimum step size for X(";I;")";
190 INPUT S0(I)
200 Next I
210 ! Calculate initial function value
220 ! Store as BestF value in var B
230 CALL MyFx(X1, N, B)
240 ! Calculate LastBestF L
250 IF B > 0 THEN L = B + 100 ELSE L = B - 100
260 I0 = 0
270 ! ********* MAIN LOOP ***********
280 REPEAT
290 I0 = I0 + 1
300 For I = 1 To N
310 X(I) = X1(I)
320 Next I
330 For I = 1 To N
340 M(I) = 0
350 LOOP
360 X0 = X1(I)
370 X1(I) = X0 + S9(I)
380 CALL MyFx(X1, N, F)
390 If F < B Then
400 B = F
410 M(I) = 1
420 Else
430 X1(I) = X0 - S9(I)
440 CALL MyFx(X1, N, F)
450 If F < B Then
460 B = F
470 M(I) = 1
480 Else
490 X1(I) = X0
500 LEAVE
510 End If
520 End If
530 END LOOP
540 Next I
550 ! moved in any direction?
560 M0 = 1
570 For I = 1 To N
580 If M(I) = 0 Then
590 M0 = 0
600 Exit I
610 End If
620 Next I
630 If M0 = 1 Then
640 For I = 1 To N
650 D(I) = X1(I) - X(I)
660 Next I
670 L8 = 1
680 CALL LSrch(X, N, L8, D, 0.1, 0.0001, K0)
690 If K0 = 1 Then
700 For I = 1 To N
710 X1(I) = X(I) + L8 * D(I)
720 Next I
730 End If
740 End If
750 CALL MyFx(X1, N, B)
760 ! reduce the step size for the dimensions that had no moves
770 For I = 1 To N
780 If M(I) = 0 Then S9(I) = S9(I) / 2
790 Next I
800 ! show results
810 DISP "For iteration ";I0
820 For I = 1 To N
830 DISP "For X(";I;") Value = ";X(I);" and step size = ";S9(I)
840 Next I
850 DISP
855 WAIT 2
860 L = B
870 S1 = 1
880 For I = 1 To N
890 If S9(I) >= S0(I) Then
900 S1 = 0
910 Exit I
920 End If
930 Next I
940 UNTIL S1=1
950 !************ END OF MAIN LOOP *************
960 DISP "******** FINAL RESULTS *********"
970 DISP "Results for optimum point"
980 For I = 1 To N
990 DISP "For X(";I;") Value = ";X(I);" and step size = ";S9(I) @ PAUSE
1000 Next I
1010 DISP "Function value = ";B @ PAUSE
1020 DISP "Number of iterations = ";I0
1030 !
1040 !
1050 SUB MyFx(X(), N, R)
1060 ! R = 100 * (X(1) ^ 2 - X(2)) ^ 2 + (1 - X(1)) ^ 2
1070 R = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2
1080 End SUB
1090 !
1100 SUB MyFxEx(N, X(), D(), L, R)
1110 Dim X0(N)
1120 For I = 1 To N
1130 X0(I) = X(I) + L * D(I)
1140 Next I
1150 CALL MyFx(X0, N, R)
1160 End SUB
1170 !
1180 Sub GetGrads(X(), N, D(), D0)
1190 D0 = 0
1200 For I = 1 To N
1210 X0 = X(I)
1220 H = 0.01
1230 If Abs(X0) > 1 Then H = H * X0
1240 X(I) = X0 + H
1250 CALL MyFx(X, N, F1)
1260 X(I) = X0 - H
1270 CALL MyFx(X, N, F2)
1280 X(I) = X0
1290 D(I) = (F1 - F2) / 2 / H
1300 D0 = D0 + D(I) ^ 2
1310 Next I
1320 D0 = Sqr(D0)
1330 End Sub
1340 !
1350 SUB LSrch(X(), N, L, D(), S9, S0, R)
1360 CALL MyFxEx(N, X, D, L, F1)
1370 REPEAT
1380 CALL MyFxEx(N, X, D, L + S9, F2)
1390 If F2 < F1 Then
1400 F1 = F2
1410 L = L + S9
1420 Else
1430 CALL MyFxEx(N, X, D, L - S9, F2)
1440 If F2 < F1 Then
1450 F1 = F2
1460 L = L - S9
1470 Else
1480 ! reduce search step size
1490 S9 = S9 / 10
1500 End If
1510 End If
1520 Until S9 < S0
1530 R = 1
1540 End SUB
The above listing uses the JPCROM to employ more structured constructs for loops and decision-making. The above listing uses the REPEAT-UNTIL and LOOP-ENDLOOP loops and also the IF-THEN-ELSE-ENDIF enhanced version of the basic IF statement. These program flow control statements limits or eliminates using the GOTO statements and make the programs easier to read and update. Here is the version of the listing that is indented and void of line numbers. I am including it here since it is easier to read than the above listing although it will not run on either the HP-71B handheld computer or the EMU71 emulator..   
! Directional Search using Hooke and Jeeves method
!
! JPCROM-BASED VERSION 1.0
! USES JPC ROM TO SUPPORT IF-ELSE-END IF AND LOOP-END LOOP STATEMENTS
!
!---------------------------- MAIN PROGRAM --------------------
!
N = 2 ! number of variables
DIM X(N), X1(N), M(N)
DIM S9(N), S0(N), D(N)
!
DISP "Hooke-Jeeves Directional Search Method"
For I = 1 To N
  DISP "Enter value for X(";I;")";
  INPUT X1(I)
  DISP "Enter step size for X(";I;")";
  INPUT S9(I)
  DISP "Enter minimum step size for X(";I;")";
  INPUT S0(I)
Next I  
! Calculate initial function value
! Store as BestF value in var B
CALL MyFx(X1, N, B)
! Calculate LastBestF L
IF B > 0 THEN L = B + 100 ELSE L = B - 100
I0 = 0 
! ********* MAIN LOOP ***********
REPEAT
  I0 = I0 + 1  
  For I = 1 To N
    X(I) = X1(I)
  Next I  
  For I = 1 To N
    M(I) = 0
    LOOP
      X0 = X1(I)
      X1(I) = X0 + S9(I)
      CALL MyFx(X1, N, F)
      If F < B Then
        B = F
        M(I) = 1
      Else
        X1(I) = X0 - S9(I)
        CALL MyFx(X1, N, F)
        If F < B Then
          B = F
          M(I) = 1
        Else
          X1(I) = X0
          LEAVE
        End If
      End If
    END LOOP
  Next I  
  ! moved in any direction?
  M0 = 1
  For I = 1 To N
    If M(I) = 0 Then
      M0 = 0
      Exit I
    End If
  Next I  
  If M0 = 1 Then
    For I = 1 To N
      D(I) = X1(I) - X(I)
    Next I
    L8 = 1
    CALL LSrch(X, N, L8, D, 0.1, 0.0001, K0)
    If K0 = 1 Then
      For I = 1 To N
        X1(I) = X(I) + L8 * D(I)
      Next I
    End If
  End If
  CALL MyFx(X1, N, B)
  ! reduce the step size for the dimensions that had no moves
  For I = 1 To N
    If M(I) = 0 Then S9(I) = S9(I) / 2
  Next I
  ! show results
  DISP "For iteration ";I0
  For I = 1 To N
    DISP "For X(";I;") Value = ";X(I);" and step size = ";S9(I)
  Next I
  DISP
  WAIT 2
  L = B  
  S1 = 1
  For I = 1 To N
    If S9(I) >= S0(I) Then
      S1 = 0
      Exit I
    End If
  Next I  
UNTIL S1=1
!************ END OF MAIN LOOP *************
DISP "******** FINAL RESULTS *********"
DISP "Results for optimum point"
For I = 1 To N
  DISP "For X(";I;") Value = ";X(I);" and step size = ";S9(I) @ PAUSE
Next I
DISP "Function value = ";B @ PAUSE
DISP "Number of iterations = ";I0
!
!
SUB MyFx(X(), N, R)
  ! R = 100 * (X(1) ^ 2 - X(2)) ^ 2 + (1 - X(1)) ^ 2
  R = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2
End SUB
!
SUB MyFxEx(N, X(), D(), L, R)
  Dim X0(N) 
  For I = 1 To N
    X0(I) = X(I) + L * D(I)
  Next I
  CALL MyFx(X0, N, R)
End SUB
!
Sub GetGrads(X(), N, D(), D0)
  D0 = 0
  For I = 1 To N
    X0 = X(I)
    H = 0.01
    If Abs(X0) > 1 Then H = H * X0
    X(I) = X0 + H
    CALL MyFx(X, N, F1)
    X(I) = X0 - H
    CALL MyFx(X, N, F2)
    X(I) = X0
    D(I) = (F1 - F2) / 2 / H
    D0 = D0 + D(I) ^ 2
  Next I
  D0 = Sqr(D0)
End Sub
!
SUB LSrch(X(), N, L, D(), S9, S0, R)  
  CALL MyFxEx(N, X, D, L, F1)  
  REPEAT
    CALL MyFxEx(N, X, D, L + S9, F2)
    If F2 < F1 Then
      F1 = F2
      L = L + S9
    Else
      CALL MyFxEx(N, X, D, L - S9, F2)
      If F2 < F1 Then
        F1 = F2
        L = L - S9
      Else
        ! reduce search step size
        S9 = S9 / 10
      End If
    End If
  Until S9 < S0
  R = 1
End SUB

Here is a version of the first listing that DOES NOT rely on the JPCROM. This version should run on the basic HP-71B with no Math or JPCROM modules.

10 ! Directional Search using Hooke and Jeeves method
50 !
60 !---------------------------- MAIN PROGRAM --------------------
70 !
80 N = 2 ! number of variables
90 DIM X(N), X1(N), M(N)
100 DIM S9(N), S0(N), D(N)
110 !
120 DISP "Hooke-Jeeves Directional Search Method"
130 For I = 1 To N
140 DISP "Enter value for X(";I;")";
150 INPUT X1(I)
160 DISP "Enter step size for X(";I;")";
170 INPUT S9(I)
180 DISP "Enter minimum step size for X(";I;")";
190 INPUT S0(I)
200 Next I
210 ! Calculate initial function value
220 ! Store as BestF value in var B
230 CALL MyFx(X1, N, B)
240 ! Calculate LastBestF L
250 IF B > 0 THEN L = B + 100 ELSE L = B - 100
260 I0 = 0
270 ! ********* MAIN LOOP ***********
280 ! REPEAT
290 I0 = I0 + 1
300 For I = 1 To N
310 X(I) = X1(I)
320 Next I
330 For I = 1 To N
340 M(I) = 0
350 ! LOOP
360 X0 = X1(I)
370 X1(I) = X0 + S9(I)
380 CALL MyFx(X1, N, F)
390 If F >= B Then 420
400 B = F
410 M(I) = 1
415 GOTO 520
420 ! Else
430 X1(I) = X0 - S9(I)
440 CALL MyFx(X1, N, F)
450 If F >= B Then 490
460 B = F
470 M(I) = 1
475 GOTO 510
480 ! Else
490 X1(I) = X0
500 GOTO 540 ! LEAVE
510 !End If
520 !End If
530 GOTO 350 ! END LOOP
540 Next I
550 ! moved in any direction?
560 M0 = 1
570 For I = 1 To N
580 If M(I) <> 0 Then 610
590 M0 = 0
600 I = N
610 !End If
620 Next I
630 If M0 <> 1 Then 740
640 For I = 1 To N
650 D(I) = X1(I) - X(I)
660 Next I
670 L8 = 1
680 CALL LSrch(X, N, L8, D, 0.1, 0.0001, K0)
690 If K0 <> 1 Then 730
700 For I = 1 To N
710 X1(I) = X(I) + L8 * D(I)
720 Next I
730 !End If
740 !End If
750 CALL MyFx(X1, N, B)
760 ! reduce the step size for the dimensions that had no moves
770 For I = 1 To N
780 If M(I) = 0 Then S9(I) = S9(I) / 2
790 Next I
800 ! show results
810 DISP "For iteration ";I0
820 For I = 1 To N
830 DISP "For X(";I;") Value = ";X(I);" and step size = ";S9(I)
840 Next I
850 DISP
855 WAIT 2
860 L = B
870 S1 = 1
880 For I = 1 To N
890 If S9(I) < S0(I) Then 920
900 S1 = 0
910 I = N
920 !End If
930 Next I
940 IF S1=0 THEN 280
950 !************ END OF MAIN LOOP *************
960 DISP "******** FINAL RESULTS *********"
970 DISP "Results for optimum point"
980 For I = 1 To N
990 DISP "For X(";I;") Value = ";X(I);" and step size = ";S9(I) @ PAUSE
1000 Next I
1010 DISP "Function value = ";B @ PAUSE
1020 DISP "Number of iterations = ";I0
1030 !
1040 !
1050 SUB MyFx(X(), N, R)
1060 ! R = 100 * (X(1) ^ 2 - X(2)) ^ 2 + (1 - X(1)) ^ 2
1070 R = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2
1080 End SUB
1090 !
1100 SUB MyFxEx(N, X(), D(), L, R)
1110 Dim X0(N)
1120 For I = 1 To N
1130 X0(I) = X(I) + L * D(I)
1140 Next I
1150 CALL MyFx(X0, N, R)
1160 End SUB
1170 !
1180 Sub GetGrads(X(), N, D(), D0)
1190 D0 = 0
1200 For I = 1 To N
1210 X0 = X(I)
1220 H = 0.01
1230 If Abs(X0) > 1 Then H = H * X0
1240 X(I) = X0 + H
1250 CALL MyFx(X, N, F1)
1260 X(I) = X0 - H
1270 CALL MyFx(X, N, F2)
1280 X(I) = X0
1290 D(I) = (F1 - F2) / 2 / H
1300 D0 = D0 + D(I) ^ 2
1310 Next I
1320 D0 = Sqr(D0)
1330 End Sub
1340 !
1350 SUB LSrch(X(), N, L, D(), S9, S0, R)
1360 CALL MyFxEx(N, X, D, L, F1)
1370 ! REPEAT
1380 CALL MyFxEx(N, X, D, L + S9, F2)
1390 If F2 >= F1 Then 1420
1400 F1 = F2
1410 L = L + S9
1415 GOTO 1510
1420 !Else
1430 CALL MyFxEx(N, X, D, L - S9, F2)
1440 If F2 >= F1 Then 1470
1450 F1 = F2
1460 L = L - S9
1465 GOTO 1500
1470 !Else
1480 ! reduce search step size
1490 S9 = S9 / 10
1500 !End If
1510 !End If
1520 IF S9 >= S0 THEN 1370
1530 R = 1
1540 End SUB

BACK

Copyright (c) Namir Shammas. All rights reserved.