Trailing-Edge
-
PDP-10 Archives
-
decuslib10-01
-
43,50110/linqd1.eng
There are 2 other files named linqd1.eng in the archive. Click here to see a list.
100' NAME--LINQD1
110'
120' DESCRIPTION--LINEAR OR QUADRATIC PROGRAMMING ALGORITHM
130'
140' SOURCE--M.P. JUHA AND A.O. CONVERSE, JUNE 1968, THAYER SCHOOL
150' OF ENGINEERING
160'
170' INSTRUCTIONS--INSERT DATA AS INDICATED BELOW BETWEEN LINES 3160 AND 3370
180'
190'
200' 1) PROBLEM TYPE: MAX OR MIN
210' 2) NUMBER OF TERMS IN OBJECTIVE FUNCTION
220' E.G., THE EXPRESSION: X1*X2 - 2*X1 - X2
230' CONSISTS OF THREE TERMS
240' 3) NUMBER OF ALLOCATION VARIABLES
250' E.G., THE EXPRESSION: X1*X2 - 2*X1 - X2
260' CONTAINS TWO ALLOCATION VARIABLES: X1 AND X2
270' 4) NUMBER OF CONSTRAINT RELATIONS ON PROBLEM
280' 5) THE TERMS OF THE OBJECTIVE FUNCTION EXPRESSED IN THE
290' FOLLOWING FORMAT:
300' [ENTRY 1],[ENTRY 2],[COEFFICIENT]
310' WHERE: ENTRY 1 IS THE NUMBER OF AN ALLOCATION VARIABLE
320' AND ENTRY 2 IS THE NUMBER OF AN ALLOCATION VARIABLE
330' E.G., THE TERM: X1*X2 WOULD BE EXPRESSED AS
340' 1,2,1
350' AND THE TERM: - 2*X1 WOULD BE 1,0,-2 OR 0,1,-2
360' [THE ZERO SIGNIFIES THAT THE TERM IS LINEAR, THAT IS
370' X0=1 BY DEFINITION]
380' 6) THE COEFFICIENTS OF THE CONSTRAINT RELATIONS WRITTEN
390' AS FOLLOWS: FOR BOTH MINIMIZATION AND MAXIMIZATION
400' PROBLEMS THE CONSTRAINTS MUST ALL BE OF THE FORM:
410' H[X1,X2,...,XN]>=0
420' AS AN EXAMPLE, FOR A MINIMIZATION PROBLEM WHERE ONE
430' OF THE CONSTRAINT RELATIONS IS
440' - 2*X1 - 3*X2 + 6 >= 0
450' THE INPUT DATA WOULD CONSIST OF THE COEFFICIENTS OF
460' X1 THROUGH XN [N=2 FOR THIS CASE] FOLLOWED BY THE
470' CONSTANT TERM. E.G., THE INPUT DATA WOULD BE:
480' -2, -3, 6
490' NOTE FOR NOMENCLATURE CONTINUE LISTING THROUGH LINE 760.
500' AN EXAMPLE IS PROVIDED FOR THE USER STARTING IN LINE 3020.
510'
520'
530' * * * * * * * MAIN PROGRAM * * * * * * * * * * *
540'
550 DIM A(51,51),M(50),Y(50)
560' *****N0MENCLATURE*****
570' T$= PROBLEM TYPE
580' W = NUMBER OF TERMS IN OBJECTIVE FUNCTION
590' Q = NUMBER OF ALLOCATION VARIABLES X(I)
600' U = NUMBER OF CONSTRAINT RELATIONS
610' R = NUMBER OF ROWS IN TABLEAU
620' C = NUMBER OF COLUMNS IN TABLEAU
630' A = WORKING MATRIX
640' M = VECTOR OF VALUES OF LAGRANGE MULTIPLIERS
650' Y = VECTOR OF VALUES OF ALLOCATION VARIABLES AND CONSTRAINTS
660' Z9 = A FLAG. IT IS SET > 1 BY THE USER WHEN HE WISHES TO STOP
670' MANUAL PIVOTING AND PRINT FINAL RESULTS
680' Z8 = A FLAG. IT IS SET > 1 BY THE USER WHEN HE WISHES TO
690' SPECIFY THE PIVOT ELEMENT MANUALLY.
700' Z7 = A FLAG. IT IS SET = 5 BY THE PROGRAM WHEN PIVOTING ON THE
710' TRANSPOSE OF THE ORIGINAL PIVOT ELEMENT.
720' Z6 = A FLAG. IT IS SET = 1 BY THE PROGRAM WHEN THE PROBLEM IS
730' LINEAR.
740' Z5 = A FLAG. IT IS SET > 1 BY THE USER WHEN PRINTOUT OF THE
750' INTERMEDIATE TABLEAUS IS DESIRED.
760'
770' *******DATA INPUT*******
780 LET Z7 = 0
790 READ T$,W,Q,U
800 LET Z6 = 0
810 PRINT"OBJECTIVE IS TO ";T$;"IMIZE THE SUM OF THE FOLLOWING TERMS:"
820 PRINT " "
830 LET R=C=U+Q+1'COMPUTE SIZE OF TABLEAU
840 MAT A = ZER(R,C)'CREATE AN EMPTY TABLEAU
850 LET F=-1 'SET FLAG TO INDICATE QUADRATIC PROGRAMMING
860'
870 FOR I = 1 TO Q
880 LET V$(I) = "X"
890 NEXT I
900 FOR I = Q+1 TO Q+U
910 LET V$(I) = "MU"
920 NEXT I
930 LET V$(Q+U+1) = "1"
940 FOR K=1 TO W 'READ IN OBJECTIVE FUNCTION TERM-BY-TERM
950 READ I,J,A(I,J)
960 IF I=0 THEN 1060 'IS TERM LINEAR? YES, GO TO 400
970 IF J=0 THEN 1060 'IS TERM LINEAR? YES, GO TO 400
980 PRINT " ";A(I,J);"* X(";I;") * X(";J;")"
990 IF I=J THEN 1020 'CHECK THAT TERM IS A SQUARE TERM
1000 LET A(J,I) = A(I,J)
1010 GO TO 1030' THE CURRENT TERM IS NOT IN THE DIAGONAL
1020 LET A(J,I)=A(I,J)+A(J,I) 'ASSURE THAT DIAGONAL TERMS ARE TIMES 2
1030 LET F=0 'INDICATE THAT A QUADRATIC TERM HAS BEEN FOUND
1040 GO TO 1080 'GO ON TO NEXT TERM IN OBJECTIVE FUNCTION
1050'READ IN COEFFICIENT OF TERM LINEAR IN X
1060 LET A(I+J,C)=A(R,I+J)=A(I,J)
1070 PRINT " ";A(I,J);"* X(";I+J;")"
1080 NEXT K 'READ IN NEXT TERM OF OBJECTIVE FUNCTION
1090 PRINT " "
1100 IF T$ = "MIN" THEN 1180
1110 FOR I9 = 1 TO Q
1120 FOR I8 = 1 TO Q+U+1
1130 LET A(I9,I8) = -A(I9,I8)
1140 NEXT I8
1150 NEXT I9
1160'
1170'
1180 PRINT "THE CONSTRAINTS ARE:"
1190'READ IN CONSTRAINT RELATIONS
1200 FOR I=1 TO U
1210 PRINT " "
1220 PRINT " H";I;" =:"
1230 FOR J=1 TO Q
1240 READ A(Q+I,J) 'READ COEFFICIENT TO X(J)
1250 PRINT " ";A(Q+I,J);"* X(";J;")"
1260 LET A(J,Q+I)=-A(Q+I,J) 'CREATE TRANSPOSE TERM
1270 NEXT J
1280 READ A(Q+I,C) 'READ CONSTANT TERM IN CONSTRAINT
1290 PRINT " ";A(Q+I,C)
1300 IF T$ = "MIN" THEN 1330
1310 LET A(R,Q+I) = A(Q+I,C)
1320 GO TO 1340
1330 LET A(R,Q+I)=-A(Q+I,C)
1340 PRINT " >= 0"
1350 NEXT I
1360 PRINT " "
1370 PRINT " "
1380 PRINT " "
1390'
1400 IF F = 0 THEN 1440' SKIP TO 660 IF THIS IS A Q.P. PROBLEM
1410 PRINT "LINEAR PROGRAMMING"
1420 PRINT" "
1430 LET Z6 = 1' SIGNIFIES A LINEAR PROGRAM
1440 PRINT"DO YOU WISH TO SELECT THE PIVOT ELEMENTS MANUALLY,YES OR NO"
1450 INPUT T1$
1460 IF T1$ = "NO" THEN 1490
1470 LET Z8 = 2
1480 GO TO 1520
1490 PRINT"DO YOU WISH PRINTOUT OF INTERMEDIATE TABLEAUS,YES OR NO"
1500 INPUT T2$
1510 IF T2$ = "NO" THEN 1530
1520 LET Z5 = 2
1530'
1540'END DATA READ IN PROCEDURE
1550'
1560'INITIALIZE EXTREMUM STORAGE MATRIX
1570 IF T$="MAX" THEN 1610
1580'PROBLEM IS MIN
1590 LET B(R,C)= 1 E 27'SET OBJECTIVE COMPARISON VERY LARGE
1600 GO TO 1660'CONTINUE ON TO SET UP EXCHANGE INDICES
1610'PROBLEM IS MAX
1620 LET B(R,C)=-1 E 27'SET OBJECTIVE COMPARISON VERY SMALL
1630'
1640'SET EXCHANGE VARIABLE INDICES, A(0,I)
1650'INDICES STORE A RECORD OF WHICH VARIABLE IS EXCHANGED
1660 FOR I = 1 TO Q
1670 LET A (0,I) = 100+I' IDENTIFIES THE ALLOCATION VARIABLES
1680 NEXT I
1690 FOR I = Q+1 TO Q+U
1700 LET A(0,I) = 200 + I' IDENTIFIES L.M. ASSOCIATED WITH CONSTRAINTS
1710 NEXT I
1720 FOR J = 1 TO R-1
1730 LET A(J,0)=J
1740 NEXT J
1750'
1760 IF Z5 < 1 THEN 1830' NO INTERMEDIATE PRINTOUT
1770 GO SUB 2090' TO PRINT EXPLANATION OF INITIAL TABLEAU
1780 MAT PRINT V$,
1790 IF Z5 < 1 THEN 1830' NO INTERMEDIATE PRINTOUT
1800'
1810 MAT PRINT A
1820 IF Z8 > 1 THEN 2170' TO SELECT PIVOT MANUALLY
1830 LET P = 1E27' SET FLAG AND SEARCH FOR PIVOT
1840 FOR I = 1 TO R-1
1850 IF A(I,C) >= 0 THEN 1960
1860 IF Z6 = 1 THEN 1880
1870 GO SUB 1990' Q.P. SEARCH FOR PIVOT
1880 FOR J = 1 TO C-1
1890 IF A(I,J) <=0 THEN 1950
1900 LET X = A(R,J)/A(I,J)
1910 IF X >=P THEN 1950
1920 LET P = X
1930 LET K = I
1940 LET L = J
1950 NEXT J
1960 NEXT I
1970 IF P < 1E27 THEN 2210
1980 GO TO 2650' NO PIVOT
1990' *** Q.P. SEARCH FOR PIVOT SUBROUTINE ***
2000'
2010 FOR I = 1 TO R-1
2020 IF A(I,0) > 100 THEN 2060' HAS VARIABLE BEEN EXCHANGED
2030 IF A(I,I) =0 THEN 2060' IS DIAGONAL ELEMENT EMPTY
2040 LET K = L = I' NOTE PIVOT LOCATION
2050 GO TO 2210' GO TO EXCHANGE ALGORITHM
2060 NEXT I
2070 RETURN
2080'
2090' PRINT SUBROUTINE
2100 PRINT" IN THE FIRST TABLEAU, THE FIRST";Q;" COLUMNS ARE ASSOCIATED"
2110 PRINT" WITH THE CORRESPONDING";Q;"ALLOCATION VARIABLES. THE NEXT"
2120 PRINT U;"COLUMNS ARE ASSOCIATED WITH THE ";U;"LAGRANGE MULTTIPLIERS"
2130 PRINT" ASSOCIATED WITH THE CORRESPONDING CONSTRAINTS. "
2140 RETURN
2150'
2160 IF Z8 < 1 THEN 2210' USER WILL NOT CHOOSE PIVOT
2170 PRINT " INPUT PIVOT LOCATION AND Z9 "
2180 INPUT K,L,Z9
2190 IF Z9 > 1 THEN 2650' USER HAS DECIDED TO STOP AND PRINT RESULTS
2200'
2210' ***** EXCHANGE ALGORITHM *****
2220 PRINT "PIVOT LOCATION";K;L
2230 LET X=A(0,L)
2240 LET A(0,L)=A(K,0)
2250 LET A(K,0)=X
2260 LET X=A(K,L)
2270 REM EXCHANGE ALL ENTRIES NOT IN PIVOT ROW OR COLUMN
2280 FOR I=1 TO R
2290 IF I=K THEN 2340
2300 FOR J=1 TO C
2310 IF J=L THEN 2330
2320 LET A(I,J)=A(I,J)-A(I,L)*A(K,J)/X
2330 NEXT J
2340 NEXT I
2350 REM CALCULATE PIVOT COLUMN VALUES
2360 FOR I=1 TO R
2370 LET A(I,L)=A(I,L)/X
2380 NEXT I
2390 REM CALCULATE PIVOT ROW VALUES
2400 FOR J=1 TO C
2410 LET A(K,J)=-A(K,J)/X
2420 NEXT J
2430 REM CALCULATE PIVOT VALUE
2440 LET A(K,L)=1/X
2450'FINISHED EXCHANGING VARIABLES
2460 IF K = L THEN 2540' SKIP TO 1625 IF PIVOT IS IN DIAGONAL
2470 IF Z7 = 5 THEN 2540' GO TO 1625 IF TRANSPOSED PIVOT HAS BEEN DONE
2480 LET Z7 = K' TEMPORARY STORAGE
2490 LET K = L' TRANSPOSING
2500 LET L = Z7
2510 LET Z7 = 5' FLAG TO SHOW WHEN TRANSPOSED PIVOT HAS BEEN DONE
2520 GO TO 2210' TO PERFORM TRANSPOSE PIVOT
2530'
2540 LET Z7 = 0' FLAG
2550 GO TO 1790
2560'
2570'
2580'PROBLEM END
2590'
2600'OUTPUT RESULTS
2610'
2620'DETERMINE VALUES OF ALLOCATION VARIABLES, LAGRANGE MULTIPLIERS
2630' AND CONSTRAINT RELATIONS.
2640'
2650 MAT M=ZER(Q+U) 'SET UP LAGRANGE MULTIPLIER VECTOR
2660 MAT Y=ZER(Q+U) 'SET UP ALLOCATION VARIABLE/CONSTRAINT VECTOR
2670'
2680 FOR I=1 TO R-1
2690 IF A(I,0) > 200 THEN 2780
2700 IF A(I,0)>100 THEN 2760
2710 IF A(I,0) > Q THEN 2740
2720 LET M(A(I,0)) = A(I,C)
2730 GO TO 2850 'GO ON TO NEXT ROW IN TABLEAU
2740 LET Y(A(I,0)) = A(I,C)
2750 GO TO 2850
2760 LET Y(A(I,0) -100) = A(I,C)
2770 GO TO 2850 'GO ON TO NEXT ROW IN TABLEAU
2780 LET M(A(I,0)-200) = A(I,C)
2790 GO TO 2850
2800'ROW OF TABLEAU IS IN CONSTRAINT RELATION PARTITION
2810 IF A(I,0)<100 THEN 2840
2820 LET M(I)=A(I,C) 'LAGRANGE MULTIPLIER, SAVE IT
2830 GO TO 2850 'GO ON TO NEXT ROW IN TABLEAU
2840 LET Y(I)=A(I,C) 'SAVE CONSTRAINT VALUE
2850 NEXT I 'GO TO NEXT ROW IN TABLEAU
2860'
2870 PRINT"ALLOCATION VARIABLES AND ASSOCIATED LAGRANGE MULTIPLIERS"
2880 PRINT " "
2890 FOR I=1 TO Q
2900 PRINT "X";I;"=";Y(I),"MU";I;"=";M(I)
2910 NEXT I
2920 PRINT " "
2930 PRINT "CONSTRAINT RELATIONS AND ASSOCIATED LAGRANGE MULTIPLIERS"
2940 PRINT " "
2950 FOR I=Q+1 TO Q+U
2960 PRINT "H";I-Q;"=";Y(I),"MU";I;"=";M(I)
2970 NEXT I
2980'
2990 PRINT " "
3000 PRINT "OBJECTIVE FUNCTION =";A(R,C)/2
3010'
3020'*****EXAMPLE PROBLEM*****
3030'
3040' MINIMIZE: 2*X1^2 + 2*X2^2 -X1*X2 -10*X1 -20*X2
3050'
3060' BY SELECTION OF: X1, X2
3070'
3080' SUBJECT TO:
3090'
3100' H1 = 8 - X2 >= 0
3110'
3120' H2 = 10 - X1 - X2 >= 0
3130'
3140' INPUT DATA FOR ALGORITHM
3150'
3160'PROBLEM TYPE
3170 DATA MIN
3180'NUMBER OF TERMS IN OBJECTIVE FUNCTION
3190 DATA 5
3200'NUMBER OF ALLOCATION VARIABLES
3210 DATA 2
3220'NUMBER OF CONSTRAINT RELATIONS
3230 DATA 2
3240'COEFFICIENTS OF OBJECTIVE FUNCTION---TERM-BY-TERM
3250' FIRST TERM
3260 DATA 1,1,2
3270' SECOND TERM
3280 DATA 2,2,2
3290' THIRD TERM
3300 DATA 1,2,-1
3310 DATA 1,0,-10
3320 DATA 2,0,-20
3330' FIRST CONSTRAINT
3340 DATA 0,-1,8
3350' SECOND CONSTRAINT
3360 DATA -1,-1,10
3370'*****END OF DATA *****
3380 END