hp41programs

Discrete Optimization

Discrete Constrained Optimization for the HP-41


Overview
 

 1°)  2 Dimensional Problems
 2°)  3 Dimensional Problems
 3°)  n Dimensional Problems  ( X-Functions Module required )
 

-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 )
 

Notes:

-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...