Discrete Constrained Optimization for the HP-41
Overview
1°) 2 Dimensional Problems
2°) 3 Dimensional Problems
3°) n Dimensional Problems ( X-Functions
Module required )
4°) Simplification ( n < 10 )
Latest Update: paragraph 4°)
-The following programs search the integers x1 ,
....... , xn that maximize a function f ( x1
, ....... , xn ) ( Seek the maximum of -f to find the minimum
of f )
provided constraints of the type x1 >=
x10 , ....... , xn >= xn0
( or x1 <= x10 , ....... ,
xn <= xn0 )
and I (
x1 , ....... , xn )
are satisfied
where I ( x1 , ....... , xn ) is one or several inequalities.
-These constraints must define a domain D and all the points ( x1 , ....... , xn ) of this domain are successively tested.
-For example, in dimension 2, we start with x = x0 and
y = y0
-Then, x is replaced by x+h ( h is an increment that you specify
- usually +1 or -1 )
-This increment is added to x until ( x, y ) falls outside D. x is then
replaced by x0 and h is added to y.
-The process is repeated until ( x , y ) is again outside D at this step.
-Therefore, none of these routines should be used blindly !
-Here are a few examples for 2 dimensional problems to visualize when it
works and when it doesn't work:
y
*
*
*
*
*
YES!
*
*
*
*
* * * *
* * * x
y
*
*
*
*
*
*
*
*
*
YES!
*
*
*
*
*
*
* * * *
* x
y
*
*
*
* *
*
*
*
*
NO! but it will work if you swap x and y
*
*
*
*
* * *
* * * * *
* * * * * x
y
*
* *
A *
* B
*
*
*
*
NO! the points above the segment [AB] would be never tested.
*
*
*
*
* * *
* * x
1°) 2 Dimensional Problems
Data Registers: • R00 = f name ( Registers R00 thru R03 are to be initialized before executing "MAX2" )
• R01 = x0 , x
R04 = max(f)
R06: temp
R07 = xSol1 ...............
• R02 = y0 , y
R05 = bbb.eee = control number of the solution(s)
R08 = ySol1 ...............
• R03 = h = increment ( usually +1 or -1 )
Flag: F01
Subroutine: A program that takes x &
y in registers R01 & R02 , sets flag F01 if ( x , y ) is outside the
domain D
and returns f(x,y) in X-register
| 01 LBL "MAX2" 02 RCL 01 03 STO 06 04 E99 05 CHS 06 STO 04 07 7 08 STO 05 09 CF 01 10 GTO 03 11 LBL 01 12 X<>Y 13 STO 04 |
14 RCL 01 15 STO 07 16 RCL 02 17 STO 08 18 7 19 STO 05 20 LBL 02 21 RCL 03 22 ST+ 01 23 LBL 03 24 XEQ IND 00 25 FS?C 01 26 GTO 04 |
27 RCL 04 28 X>Y? 29 GTO 02 30 X#Y? 31 GTO 01 32 3 33 ST+ 05 34 RCL 02 35 STO IND 05 36 DSE 05 37 RCL 01 38 STO IND 05 39 GTO 02 |
40 LBL 04 41 RCL 01 42 RCL 06 43 X=Y? 44 GTO 05 45 STO 01 46 RCL 03 47 ST+ 02 48 GTO 03 49 LBL 05 50 1 51 RCL 05 52 + |
53 E3 54 / 55 7 56 + 57 STO 05 58 RCL 04 59 RCL 08 60 RCL 07 61 END |
( 89 bytes / SIZE maximum - at least 009 )
| STACK | INPUTS | OUTPUTS |
| T | / | bbb.eee |
| Z | / | max(f) |
| Y | / | y1 |
| X | / | x1 |
where bbb.eee = control number of the solution(s) = R05
Example1: Find the non-negative integers x , y such that f(x,y) = 7 x + 5 y is maximum with the constraints:
7 x + 10 y <=
70
8 x + 7 y
<= 56
11 x + 8 y <=
66
-Load the routine:
LBL "T" RCL 02 70
8
*
SF 01 RCL 02
66
7
*
RCL 01 10
X<Y? *
+
RCL 01 6
X<Y?
*
+
7
*
SF 01 RCL 02
56
11
*
SF 01 RCL
02 RTN
*
+
RCL 01 7
X<Y? *
+
RCL 01 5
"T" ASTO 00
0 STO 01 STO 02
1 STO 03
XEQ "MAX2" >>>>
x = 4
---Execution time = 86s---
RDN y = 3
RDN max(f) = 43
RDN 7.008
-So, there is a unique solution ( x , y ) = ( 4 , 3 ) and
max(f) = 43
Example2: Same constraints with f(x,y) = 3 x + 2 y
-The subroutine becomes:
LBL "T" RCL 02 70
8
*
SF 01 RCL 02
66
3
+
RCL 01 10
X<Y? *
+
RCL 01 6
X<Y?
*
RTN
7
*
SF 01 RCL 02
56
11
*
SF 01 RCL
02
*
+
RCL 01 7
X<Y? *
+
RCL 01 ST+ X
0 STO 01 STO 02
XEQ "MAX2" >>>>
x = 6
---Execution time = 82s---
RDN y = 0
RDN max(f) = 18
RDN 7.010
We have R07 = 6 , R08 = 0 and R09 = 4
, R10 = 3
-Therefore, the problem has 2 solutions: ( 6 , 0 ) and ( 4 , 3 ).
They give max( 3 x + 2 y ) = 18
Example3:
3 x + y >= 12
x <= 10 Find the
integers x , y that minimize f(x,y)
= 5 x + 4 y
x + y >= 8
y <= 12
3 x + 5 y >= 30
-We must maximize -f(x,y) = - 5 x - 4 y
LBL "T" RCL 02
SF 01 8
3
*
SF 01 RCL
02 CHS
RCL 01 +
RCL 01 X>Y?
*
+
RCL 01 4
RTN
3
12
RCL 02 SF 01
RCL 02 30
5
*
*
X>Y? +
RCL 01 5
X>Y? *
+
"T" ASTO 00
10 STO 01
12 STO 02
-1 STO 03
( we must choose h < 0 )
XEQ "MAX2" >>>>
x = 2
---Execution time = 2m57s---
RDN y = 6
RDN max(-f) = -34
RDN 7.008
-Here, min(f) = -max(-f) = 34 is reached for a unique point:
( x , y ) = ( 2 , 6 )
Example4: Constraints: x <= - y4 + 9 y3 - 21 y2 - y + 61 with x , y non-negative integers. We want f(x,y) = x y maximum
LBL "T" -
-
-
+
RCL 02
RCL 01 RCL 02
RCL 02 RCL 02
X<Y? *
9
*
*
*
SF 01 RTN
RCL 02 21
1
61
RCL 01
"T" ASTO 00
1 STO 01 STO 02
Since we seek a maximum product x y, x & y are surely > 0
1 STO 03
XEQ "MAX2" >>>>
x = 41
---Execution time = 3m56s---
RDN y =
4
RDN max(f) = 164
RDN 7.008
-The unique solution is ( 41 , 4 ) and fmax = 164
2°) 3 Dimensional Problems
Data Registers: • R00 = f name ( Registers R00 thru R04 are to be initialized before executing "MAX3" )
• R01 = x0 , x
R05 = max(f)
R07 , R08: temp
R09 = xSol1 ...............
• R02 = y0 , y
R06 = bbb.eee = control number of the solution(s)
R10 = ySol1 ...............
• R03 = z0 , z
R11 = zSol1 ...............
• R04 = h = increment ( usually +1 or -1 )
Flag: F01
Subroutine: A program that takes x ,
y , z in registers R01 , R02 , R03 , sets flag F01 if ( x , y , z ) is outside
the domain D
and returns f(x,y,z) in X-register
| 01 LBL "MAX3" 02 RCL 01 03 STO 07 04 RCL 02 05 STO 08 06 E99 07 CHS 08 STO 05 09 9 10 STO 06 11 CF 01 12 GTO 03 13 LBL 01 14 X<>Y 15 STO 05 16 RCL 01 |
17 STO 09 18 RCL 02 19 STO 10 20 RCL 03 21 STO 11 22 9 23 STO 06 24 LBL 02 25 RCL 04 26 ST+ 01 27 LBL 03 28 XEQ IND 00 29 FS?C 01 30 GTO 04 31 RCL 05 32 X>Y? |
33 GTO 02 34 X#Y? 35 GTO 01 36 5 37 ST+ 06 38 RCL 03 39 STO IND 06 40 DSE 06 41 RCL 02 42 STO IND 06 43 DSE 06 44 RCL 01 45 STO IND 06 46 GTO 02 47 LBL 04 48 RCL 01 |
49 RCL 07 50 X=Y? 51 GTO 05 52 STO 01 53 RCL 04 54 ST+ 02 55 GTO 03 56 LBL 05 57 RCL 02 58 RCL 08 59 X=Y? 60 GTO 06 61 STO 02 62 RCL 04 63 ST+ 03 64 GTO 03 |
65 LBL 06 66 RCL 06 67 2 68 + 69 E3 70 / 71 9 72 + 73 STO 06 74 SIGN 75 RCL 05 76 RCL 11 77 RCL 10 78 RCL 09 79 END |
( 112 bytes / SIZE maximum - at least 012 )
| STACK | INPUTS | OUTPUTS |
| T | / | max(f) |
| Z | / | z1 |
| Y | / | y1 |
| X | / | x1 |
| L | / | bbb.eee |
where bbb.eee = control number of the solution(s) = R06
Example: Find the non-negative
integers x , y , z such that f(x,y,z) = 2 x + 3 y + 5 z
is maximum provided x2/41 + y2/37 + z2/29
<= 1
LBL "T" 41
37
X^2 X>Y?
RCL 02 RCL 03
RTN
1
/
/
29 SF 01
3
5
RCL 01 RCL 02
+
/ RCL 01
*
*
X^2
X^2
RCL 03 +
ST+ X +
+
"T" ASTO 00
0 STO 01 STO 02
STO 03
1 STO 04 ( h = 1 )
XEQ "MAX3" >>>>
x = 3
---Execution time = 5m15s---
RDN y = 4
RDN z = 3
RDN max(f) = 33
LASTX 9.017
and we have:
R09 = 3 , R10 = 4 , R11 = 3
( this first solution was in the stack at the end )
R12 = 2 , R13 = 3 , R14 = 4
R15 = 1 , R16 = 2 , R17 = 5
-So, there are 3 solutions: ( 3 , 4 , 3 ) ( 2 , 3 ,
4 ) ( 1 , 2 , 5 ) all giving max( 2 x + 3 y + 5 z
) = 33
3°) N Dimensional Problems
Data Registers: • R00 = f name ( Registers R00 and R09 thru R10+n are to be initialized before executing "MAX" )
R01 = max(f)
R03 to R06: temp
R02 = bbb.eee = control number of the solution(s)
R07 & R08: unused
• R09 = n
• R10 = h = increment ( usually +1 or -1 )
• R11 = x10 , x1
R11+n = x10
R10+2n = x1 (sol1)
• R12 = x20 , x2
R12+n = x20
R11+2n = x2 (sol1)
.....................
.................
.............................
• R10+n = xn0 , xn
R09+2n = xn-10 R09+3n
= xn (sol1)
Flag: F01
Subroutine: A program that takes x1
, ....... , xn in registers R11 , ........... , R10+n
, sets flag F01 if ( x1 , ....... , xn ) is outside
the domain D
and returns f( x1 , ....... , xn ) in X-register
| 01 LBL "MAX" 02 E99 03 CHS 04 STO 01 05 RCL 09 06 ST+ X 07 10 08 + 09 STO 02 10 STO 06 11 RCL 09 12 E6 13 / 14 11 15 + 16 STO 05 17 RCL 09 18 LASTX 19 + |
20 E3 21 / 22 + 23 REGMOVE 24 CF 01 25 GTO 03 26 LBL 00 27 X<>Y 28 STO 01 29 RCL 06 30 STO 02 31 LBL 01 32 RCL 02 33 E3 34 / 35 RCL 05 36 + 37 REGMOVE 38 LBL 02 |
39 RCL 10 40 ST+ 11 41 LBL 03 42 XEQ IND 00 43 FS?C 01 44 GTO 04 45 RCL 01 46 X>Y? 47 GTO 02 48 X#Y? 49 GTO 00 50 RCL 09 51 ST+ 02 52 GTO 01 53 LBL 04 54 RCL 05 55 INT 56 STO 03 57 RCL 06 |
58 DSE X 59 E3 60 / 61 + 62 RCL 09 63 + 64 STO 04 65 LBL 05 66 RCL IND 03 67 RCL IND 04 68 X=Y? 69 GTO 06 70 STO IND 03 71 ISG 03 72 CLX 73 RCL 10 74 ST+ IND 03 75 GTO 03 76 LBL 06 |
77 ISG 03 78 CLX 79 ISG 04 80 GTO 05 81 RCL 01 82 RCL 02 83 RCL 09 84 + 85 DSE X 86 E3 87 / 88 RCL 06 89 + 90 STO 02 91 END |
( 133 bytes / SIZE maximum )
| STACK | INPUTS | OUTPUTS |
| Y | / | max(f) |
| X | / | bbb.eee |
where bbb.eee = control number of the solution(s) = R02
Example: Find the non-negative
integers x , y , z , t that maximize f(x,y,z,t) = x y2
z3 t4 under the constraint x2
+ y2 + z2 + t2 < 100
LBL "T" X^2
+
100
RCL 12 ST* Y
X^2
RCL 11 +
RCL 14 X<=Y?
X^2
X^2 X^2
X^2
RCL 13 X^2
SF 01 *
*
*
RCL 12 X^2
+
RCL 11 RCL 13
RCL 14 RTN
"T" ASTO 00
-We have n = 4 unknowns
4 STO 09
-The increment h = 1
1 STO 10
-We can choose x0 = y0 = z0 = t0
= 1 1 STO 11 STO 12
STO 13 STO 14
XEQ "MAX" >>>>
bbb.eee = 18.021
---Execution time = 1h02m---
( Obviously not very fast... )
X<>Y max(f) = 14406000
-We find R18 = 3 , R19 = 4 , R20 = 5 , R21 = 7
-So, the unique solution is ( x , y , z , t ) = ( 3 , 4 , 5 , 7 )
4°) Simplification ( n < 10 )
-In the following program, we assume that:
x1 , ....... , xn are non-negative integers
h = +1
-Thus, we can save a few bytes.
Data Registers: • R00 = f name ( Register R00 is to be initialized before executing "MAX" )
R01 = x1 R14 to R13+n = sol1
R02 = x2 R14+n to R13+2n = Sol2 ( if any )
.......... .......................
Rnn = xn
R10 = n R11 = max(f) R12-R13: temp
Flag: F01
Subroutine: A program that takes x1 , ....... , xn in registers R01 , ........... , Rnn , sets flag F01 if ( x1 , ....... , xn ) is outside the domain D
and returns f( x1 , ....... , xn ) in X-register
-Without an HP41CX, line 08 CLRGX may be replaced by
0 LBL 00 STO IND Y ISG Y GTO 00
| 01 LBL "MAX" 02 CF 01 03 STO 10 04 STO 13 05 E3 06 / 07 ISG X 08 CLRGX 09 E99 10 CHS 11 STO 11 12 14 |
13 STO 12 14 GTO 05 15 LBL 01 16 X<>Y 17 STO 11 18 14 19 STO 12 20 LBL 02 21 RCL 12 22 RCL 10 23 E3 24 ST/ Z |
25 X^2 26 / 27 ISG Y 28 + 29 REGMOVE 30 LBL 03 31 RCL 10 32 STO 13 33 LBL 04 34 ISG IND 13 35 LBL 05 36 XEQ IND 00 |
37 FS?C 01 38 GTO 06 39 RCL 11 40 X>Y? 41 GTO 03 42 X#Y? 43 GTO 01 44 RCL 10 45 ST+ 12 46 GTO 02 47 LBL 06 48 CLX |
49 STO IND 13 50 DSE 13 51 GTO 04 52 RCL 10 53 RCL 12 54 + 55 DSE X 56 E3 57 / 58 14 59 + 60 RCL 11 61 END |
( 96 bytes / SIZE 014+nnn+ ??? )
| STACK | INPUT | OUTPUTS |
| Y | / | bbb.eee |
| X | n | max(f) |
where bbb.eee = control number of the solution(s) & max(f) = R11
Example: Find the non-negative integers x , y , z such that f(x,y,z) = 2 x + 3 y + 5 z is maximum provided x2/41 + y2/37 + z2/29 <= 1
| 01 LBL "T" 02 CLX 03 SIGN 04 RCL 01 05 X^2 06 41 07 / 08 RCL 02 09 X^2 10 37 11 / 12 + 13 RCL 03 14 X^2 15 29 16 / 17 + 18 X>Y? 19 SF 01 20 RCL 01 21 ST+ X 22 RCL 02 23 3 24 * 25 + 26 RCL 03 27 5 28 * 29 + 30 RTN |
"T" ASTO 00 n = 3 so
3 XEQ "MAX" >>>> max(f) = 33 ---Execution time = 5m22s---X<>Y 14.022
R14 = 1 , R15 = 2 , R16 = 5
R17 = 2 , R18 = 3 , R19 = 4
R20 = 3 , R21 = 4 , R22 = 3
-So, the 3 solutions are: ( 1 , 2 , 5 ) ( 2 , 3 , 4 ) ( 3 , 4 , 3 ) all giving max( 2 x + 3 y + 5 z ) = 33
Notes:
-If we store 41 37 29 in registers R04 R05 R06
and if we replace 41 37 29 by RCL 04 RCL 05 RCL 06 in the subroutine "T",
execution time becomes 4m38s
-If, for instance, n = 3, "MAX" tests (0,0,0) (0,0,1) ..... until (0,0,z) is outside the domain D
-Then it tests (0,1,0) (0,1,1) (0,1,2) ... and so on ...
Remark:
-Always execute a SIZE maximum before executing "MAX2" "MAX3" or
"MAX"
-Even if there is eventually a unique solution, many registers may be used
during the program execution.
-Unlike the simplex method, these routines return all the solutions.
-They work for non-linear systems too.
-But execution time will limit their application to relatively "small"
sets of points...