The following program calculates the minimum point of a multi-variable function using the Powell search method.
The program prompts you to enter for each variable:
1. The function tolerance value
2. The \step tolerance value.
3. The maximum number of search cyles
4. The values for the initial set of variables
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 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. The first figure shows the input.:
The second figure shows the output:
Here is the BASIC listing:
10 ! Powell Search Optimization
20 !
30 ! USES JPCROM TO SUPPORT ADVANCED LOOPS AND
40 ! DECISION-MAKING CONSTRUCTS
50 !
60 N = 2
70 !
80 DIM X(N)
90 DIM X1(N)
100 DIM X2(N)
110 DIM S(N+1, N)
120 !
130 DISP "Powell Search Optimization"
140 INPUT "Enter function tolerance value? "; E1
150 INPUT "Enter step tolerance value? "; E2
160 INPUT "Enter maximum number of cycles? "; M9
170 For I = 1 To N
180 DISP "Enter guess for X(";I;")";
190 INPUT X(I)
200 Next I
210 !
220 J = 1
230 CALL MyFx(X, N, F1)
240 For K = 1 To N
250 X1(K) = X(K)
260 For L = 1 To N
270 S(K,L) = 0
280 IF K=L Then S(K, L) = 1
290 Next L
300 Next K
310 !
320 LOOP
330 !
340 ! reset row S(N+1,:)
350 For K = 1 To N
360 S(N + 1, K) = 0
370 Next K
380 !
390 For I = 1 To N
395 A1 = 0.5
400 CALL LSrch(X, N, A1, S, I, 0.1, 0.0001, K0)
410 For K = 1 To N
420 X(K) = X(K) + A1 * S(I, K)
430 S(N + 1, K) = S(N + 1, K) + A1 * S(I, K)
440 Next K
450 Next I
460 !
465 A1 = 0.5
470 CALL LSrch(X, N, A1, S, N + 1, 0.1, 0.0001, K0)
480 !
490 For K = 1 To N
500 X(K) = X(K) + A1 * S(N + 1, K)
510 X2(K) = X(K)
520 Next K
530 CALL MyFx(X2, N, F2)
540 !
550 DISP "F = ";F2;" ";
560 For K = 1 To N
570 DISP "X(";K;")=";X(K);" ";
580 Next K
590 DISP
600 !
610 ! test end of iterations criteria
620 If Abs(F2 - F1) < E1 Then LEAVE
630 T = 0
640 For K = 1 To N
650 T = T + (X2(K) - X1(K)) ^ 2
660 Next K
670 T = Sqr(T)
680 If T < E2 Then LEAVE
690 If J => M9 Then LEAVE
700 !
710 J = J + 1
720 !
730 ! rotate data
740 For K = 1 To N
750 X1(K) = X2(K)
760 Next K
770 F1 = F2
780 !
790 ! copy S matrix
800 For K = 1 To N
810 For L = 1 To N
820 S(K, L) = S(K + 1, L)
830 Next L
840 Next K
850 !
860 END LOOP
870 !
880 DISP "**********FINAL RESULTS************"
890 DISP "Optimum at:"
900 For I = 1 To N
910 DISP "X(";I;")=";X(I) @ PAUSE
920 Next I
930 DISP "Function value ="; F1 @ PAUSE
940 DISP "Number of iterations = ";J
950 !
960 SUB MyFx(X(), N, R)
970 ! R = 100 * (X(1) ^ 2 - X(2)) ^ 2 + (1 - X(1)) ^ 2
980 R = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2
990 End SUB
1000 !
1010 SUB MyFxEx(N, X(), S(,), I1, A1, R)
1020 DIM X0(N)
1030 For I = 1 To N
1040 X0(I) = X(I) + A1 * S(I1, I)
1050 Next I
1060 CALL MyFx(X0, N, R)
1070 End SUB
1080 !
1090 Sub GetGrad(X(), N, D5(), D6)
1100 D6 = 0
1110 For I = 1 To N
1120 X0 = X(I)
1130 H = 0.01 * (1 + ABS(X0))
1140 X(I) = X0 + H
1150 CALL MyFx(X, N, F1)
1160 X(I) = X0 - H
1170 CALL MyFx(X, N, F2)
1180 X(I) = X0
1190 D5(I) = (F1 - F2) / 2 / H
1200 D6 = D6 + D5(I) ^ 2
1210 Next I
1220 D6 = Sqr(D6)
1230 End Sub
1240 !
1250 SUB LSrch(X(), N, A1, S(,), I1, I9, I0, R)
1260 CALL MyFxEx(N, X, S, I1, A1, F1)
1270 REPEAT
1280 CALL MyFxEx(N, X, S, I1, A1 + I9, F2)
1290 If F2 < F1 Then
1300 F1 = F2
1310 A1 = A1 + I9
1320 Else
1330 CALL MyFxEx(N, X, S, I1, A1 - I9, F2)
1340 If F2 < F1 Then
1350 F1 = F2
1360 A1 = A1 - I9
1370 Else
1380 ! reduce search step size
1390 I9 = I9 / 10
1400 End If
1410 End If
1420 UNTIL I9 < I0
1430 R = 1
1440 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..
! Powell Search Optimization
!
! USES JPCROM TO SUPPORT ADVANCED LOOPS AND
! DECISION-MAKING CONSTRUCTS
!
N = 2
!
DIM X(N)
DIM X1(N)
DIM X2(N)
DIM S(N+1, N)
!
DISP "Powell Search Optimization"
INPUT "Enter function tolerance value? "; E1
INPUT "Enter step tolerance value? "; E2
INPUT "Enter maximum number of cycles? "; M9
For I = 1 To N
DISP "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
!
LOOP
!
! reset row S(N+1,:)
For K = 1 To N
S(N + 1, K) = 0
Next K
!
For I = 1 To N
A1 = 0.5
CALL LSrch(X, N, A1, S, I, 0.1, 0.0001, K0)
For K = 1 To N
X(K) = X(K) + A1 * S(I, K)
S(N + 1, K) = S(N + 1, K) + A1 * S(I, K)
Next K
Next I
!
A1 = 0.5
CALL LSrch(X, N, A1, S, N + 1, 0.1, 0.0001, K0)
!
For K = 1 To N
X(K) = X(K) + A1 * S(N + 1, K)
X2(K) = X(K)
Next K
CALL MyFx(X2, N, F2)
!
DISP "F = ";F2;" ";
For K = 1 To N
DISP "X(";K;")=";X(K);" ";
Next K
DISP
!
! test end of iterations criteria
If Abs(F2 - F1) < E1 Then LEAVE
T = 0
For K = 1 To N
T = T + (X2(K) - X1(K)) ^ 2
Next K
T = Sqr(T)
If T < E2 Then LEAVE
If J => M9 Then LEAVE
!
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
!
END LOOP
!
DISP "**********FINAL RESULTS************"
DISP "Optimum at:"
For I = 1 To N
DISP "X(";I;")=";X(I) @ PAUSE
Next I
DISP "Function value ="; F1 @ PAUSE
DISP "Number of iterations = ";J
!
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(), S(,), I1, A1, R)
DIM X0(N)
For I = 1 To N
X0(I) = X(I) + A1 * S(I1, I)
Next I
CALL MyFx(X0, N, R)
End SUB
!
Sub GetGrad(X(), N, D5(), D6)
D6 = 0
For I = 1 To N
X0 = X(I)
H = 0.01 * (1 + ABS(X0))
X(I) = X0 + H
CALL MyFx(X, N, F1)
X(I) = X0 - H
CALL MyFx(X, N, F2)
X(I) = X0
D5(I) = (F1 - F2) / 2 / H
D6 = D6 + D5(I) ^ 2
Next I
D6 = Sqr(D6)
End Sub
!
SUB LSrch(X(), N, A1, S(,), I1, I9, I0, R)
CALL MyFxEx(N, X, S, I1, A1, F1)
REPEAT
CALL MyFxEx(N, X, S, I1, A1 + I9, F2)
If F2 < F1 Then
F1 = F2
A1 = A1 + I9
Else
CALL MyFxEx(N, X, S, I1, A1 - I9, F2)
If F2 < F1 Then
F1 = F2
A1 = A1 - I9
Else
! reduce search step size
I9 = I9 / 10
End If
End If
UNTIL I9 < I0
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 ! Powell Search Optimization
20 !
60 N = 2
70 !
80 DIM X(N)
90 DIM X1(N)
100 DIM X2(N)
110 DIM S(N+1, N)
120 !
130 DISP "Powell Search Optimization"
140 INPUT "Enter function tolerance value? "; E1
150 INPUT "Enter step tolerance value? "; E2
160 INPUT "Enter maximum number of cycles? "; M9
170 For I = 1 To N
180 DISP "Enter guess for X(";I;")";
190 INPUT X(I)
200 Next I
210 !
220 J = 1
230 CALL MyFx(X, N, F1)
240 For K = 1 To N
250 X1(K) = X(K)
260 For L = 1 To N
270 S(K,L) = 0
280 IF K=L Then S(K, L) = 1
290 Next L
300 Next K
310 !
320 ! LOOP
330 !
340 ! reset row S(N+1,:)
350 For K = 1 To N
360 S(N + 1, K) = 0
370 Next K
380 !
390 For I = 1 To N
395 A1 = 0.5
400 CALL LSrch(X, N, A1, S, I, 0.1, 0.0001, K0)
410 For K = 1 To N
420 X(K) = X(K) + A1 * S(I, K)
430 S(N + 1, K) = S(N + 1, K) + A1 * S(I, K)
440 Next K
450 Next I
460 !
465 A1 = 0.5
470 CALL LSrch(X, N, A1, S, N + 1, 0.1, 0.0001, K0)
480 !
490 For K = 1 To N
500 X(K) = X(K) + A1 * S(N + 1, K)
510 X2(K) = X(K)
520 Next K
530 CALL MyFx(X2, N, F2)
540 !
550 DISP "F = ";F2;" ";
560 For K = 1 To N
570 DISP "X(";K;")=";X(K);" ";
580 Next K
590 DISP
600 !
610 ! test end of iterations criteria
620 If Abs(F2 - F1) < E1 Then GOTO 880
630 T = 0
640 For K = 1 To N
650 T = T + (X2(K) - X1(K)) ^ 2
660 Next K
670 T = Sqr(T)
680 If T < E2 Then GOTO 880
690 If J => M9 Then GOTO 880
700 !
710 J = J + 1
720 !
730 ! rotate data
740 For K = 1 To N
750 X1(K) = X2(K)
760 Next K
770 F1 = F2
780 !
790 ! copy S matrix
800 For K = 1 To N
810 For L = 1 To N
820 S(K, L) = S(K + 1, L)
830 Next L
840 Next K
850 !
860 GOTO 320 ! END LOOP
870 !
880 DISP "**********FINAL RESULTS************"
890 DISP "Optimum at:"
900 For I = 1 To N
910 DISP "X(";I;")=";X(I) @ PAUSE
920 Next I
930 DISP "Function value ="; F1 @ PAUSE
940 DISP "Number of iterations = ";J
950 !
960 SUB MyFx(X(), N, R)
970 ! R = 100 * (X(1) ^ 2 - X(2)) ^ 2 + (1 - X(1)) ^ 2
980 R = X(1) - X(2) + 2 * X(1) ^ 2 + 2 * X(1) * X(2) + X(2) ^ 2
990 End SUB
1000 !
1010 SUB MyFxEx(N, X(), S(,), I1, A1, R)
1020 DIM X0(N)
1030 For I = 1 To N
1040 X0(I) = X(I) + A1 * S(I1, I)
1050 Next I
1060 CALL MyFx(X0, N, R)
1070 End SUB
1080 !
1090 Sub GetGrad(X(), N, D5(), D6)
1100 D6 = 0
1110 For I = 1 To N
1120 X0 = X(I)
1130 H = 0.01 * (1 + ABS(X0))
1140 X(I) = X0 + H
1150 CALL MyFx(X, N, F1)
1160 X(I) = X0 - H
1170 CALL MyFx(X, N, F2)
1180 X(I) = X0
1190 D5(I) = (F1 - F2) / 2 / H
1200 D6 = D6 + D5(I) ^ 2
1210 Next I
1220 D6 = Sqr(D6)
1230 End Sub
1240 !
1250 SUB LSrch(X(), N, A1, S(,), I1, I9, I0, R)
1260 CALL MyFxEx(N, X, S, I1, A1, F1)
1270 !REPEAT
1280 CALL MyFxEx(N, X, S, I1, A1 + I9, F2)
1290 If F2 >= F1 Then 1320
1300 F1 = F2
1310 A1 = A1 + I9
1315 GOTO 1410
1320 !Else
1330 CALL MyFxEx(N, X, S, I1, A1 - I9, F2)
1340 If F2 >= F1 Then 1370
1350 F1 = F2
1360 A1 = A1 - I9
1365 GOTO 1400
1370 !Else
1380 ! reduce search step size
1390 I9 = I9 / 10
1400 !End If
1410 !End If
1420 IF I9 >= I0 THEN 1270
1430 R = 1
1440 End SUB
Copyright (c) Namir Shammas. All rights reserved.