'program simplex 'This is a QBASIC program. QBASIC is in the DOS directory. 'label start,again,fail,fail2 DIM i, infeas, j, l, lim, m, n, num, p, q, flag, flag2 AS INTEGER DIM k(1 TO 30), link(1 TO 30), bas(1 TO 30) AS INTEGER DIM r, s, denom, half, large, zero, eps, obj, opt AS DOUBLE DIM a(1 TO 30, 1 TO 30) AS DOUBLE DIM b(1 TO 30), c(1 TO 30), d(1 TO 30) AS DOUBLE DIM min(1 TO 30), x(1 TO 30) AS DOUBLE 'start executable program PRINT "This program solves the problem:" PRINT PRINT "Maximize c.x" PRINT " subject to Ax <= b, x >= 0. " PRINT PRINT "The maximum size of A is 30*30 but output format" PRINT "limitations restrict it to smaller problems if the display" PRINT "is to be interpreted easily. If mistakes are made on input " PRINT "restart by hitting CNTRL-C. Note that if an element" PRINT "of b is negative the problem is modified, adding a quantity" PRINT "of the order of 100000 to the objective function for each such" PRINT "element. The values of the variables at the optimum are unaltered." PRINT "The program is designed for data values of the order of unity. If" PRINT "numbers of magnitudes very different from unity are used, the" PRINT "results may not be reliable" PRINT PRINT 'now we start the calculation and return here for further cases start: PRINT "------------------------------------------------------------" PRINT PRINT " START NEW CASE " PRINT INPUT "type number of variables, or 0 to exit ", n IF n = 0 THEN END INPUT "type number of constraints ", m PRINT "type A matrix" FOR i = 1 TO m FOR j = 1 TO n PRINT "A("; i; ","; j; ")= "; : INPUT " ", a(i, j) NEXT j PRINT NEXT i PRINT "type b vector" FOR i = 1 TO m PRINT "b("; i; ")= "; : INPUT " ", b(i) NEXT i PRINT "type c vector" FOR i = 1 TO n PRINT "c("; i; ")= "; : INPUT " ", c(i) NEXT i large = 100000 half = -large / 2 zero = 0! eps = .00000005# 'convert to canonical form FOR i = 1 TO n bas(i) = 0 link(i) = 0 NEXT i p = n num = 1 flag2 = 0 FOR i = 1 TO m flag = 0 IF b(i) < 0! THEN flag = 1 flag2 = 1 END IF p = p + 1 bas(p) = num link(p) = 0 FOR j = 1 TO m a(j, p) = 0! NEXT j a(i, p) = 1! c(p) = 0! d(p) = 0! IF flag = 1 THEN bas(p) = 0 link(p) = p + 1 link(p + 1) = p p = p + 1 bas(p) = num FOR j = 1 TO m a(j, p) = 0! NEXT j FOR j = 1 TO p a(i, j) = -a(i, j) NEXT j b(i) = -b(i) a(i, p - 1) = -1! a(i, p) = 1! c(p) = -large d(p) = -large END IF num = num + 1! NEXT i obj = 0! FOR i = 1 TO n d(i) = c(i) NEXT i 'remove non-zero c elements for basis variables FOR i = 1 TO p IF bas(i) <> 0 THEN IF c(i) <> zero THEN FOR j = 1 TO m IF a(j, i) <> zero THEN r = c(i) / a(j, i) FOR l = 1 TO p c(l) = c(l) - r * a(j, l) NEXT l END IF NEXT j END IF END IF NEXT i PRINT PRINT "canonical form obtained by adding "; p - n; " slack variables" IF flag = 1 THEN PRINT "problem modified to handle negative right-hand sides" END IF 'the A,b and c arrays are now complete 'write out the basis again: PRINT PRINT "Sequence of variables in basis vector" FOR i = 1 TO p PRINT USING "######"; bas(i); NEXT i 'write out current tableau PRINT PRINT PRINT " A matrix and"; FOR i = 1 TO p - 2 PRINT " "; NEXT i PRINT " b vector" FOR i = 1 TO m FOR j = 1 TO p PRINT USING "####.##"; a(i, j); NEXT j PRINT USING "######.##"; b(i); PRINT NEXT i PRINT PRINT " c vector and"; FOR i = 1 TO p - 2 PRINT " "; NEXT i PRINT " -obj. fn." FOR j = 1 TO p IF ABS(c(j)) < -half THEN PRINT USING "####.##"; c(j); ELSE PRINT " ****"; END IF NEXT j PRINT USING "########.##"; obj PRINT 'check for optimality and feasibility opt = 1! FOR i = 1 TO p IF c(i) > 0! THEN opt = 0! NEXT i IF opt = 1! THEN PRINT PRINT "This is the optimum, c.x = "; PRINT USING "#########.###"; -obj; IF flag2 = 1 THEN PRINT " for the MODIFIED PROBLEM" PRINT FOR i = 1 TO p x(i) = 0! IF bas(i) <> 0 THEN j = bas(i) x(i) = b(j) / a(j, i) END IF NEXT i r = 0! FOR i = 1 TO p r = r + x(i) * d(i) NEXT i IF flag2 = 1 THEN PRINT PRINT "For unmodified problem, c.x = "; PRINT USING "#####.####"; r: END IF PRINT PRINT "The solution vector is" FOR i = 1 TO p PRINT USING "####.##"; x(i); NEXT i PRINT infeas = 0 FOR i = 1 TO p IF ((bas(i) > 0) AND (d(i) < half)) THEN IF (x(i) > eps) THEN PRINT USING "###"; i; PRINT USING "####"; bas(i) PRINT USING "####.##"; d(i) PRINT USING "##.######"; x(i) infeas = 1 j = i END IF END IF NEXT i IF infeas = 1 THEN i = j PRINT "but this is not a feasible solution, " PRINT "so the problem has NO FEASIBLE SOLUTION" END IF PRINT PRINT GOTO start END IF 'look for variables to be put into basis FOR i = 1 TO p min(i) = -1! IF (c(i) >= 0) AND (bas(i) = 0) THEN lim = 0 min(i) = large k(i) = 0 FOR j = 1 TO m FOR l = 1 TO p IF l <> j THEN IF (ABS(a(j, l)) = 1) AND (bas(l) <> 0) THEN s = a(j, l) END IF NEXT l IF (ABS(a(j, i)) > eps) THEN r = b(j) / a(j, i) IF (r > 0) OR ((r >= 0) AND (a(j, i) > 0)) OR (a(j, i) > 0) THEN lim = 1 END IF IF (r < min(i)) AND ((r > 0) OR ((r >= 0) AND (a(j, i) > 0))) THEN min(i) = r k(i) = j END IF END IF NEXT j IF link(i) <> 0 THEN l = link(i) IF bas(l) <> 0 THEN lim = 1 END IF IF (min(i) >= 0) AND (min(i) < -half) THEN PRINT "Variable "; i; " can be put into the basis, "; PRINT "benefit is "; PRINT USING "########.###"; min(i) * c(i) END IF IF lim = 0 THEN PRINT PRINT "no limit on increase of variable "; i; " SOLUTION UNBOUNDED" PRINT GOTO start END IF END IF NEXT i PRINT 'now we select the variable to be used PRINT fail: INPUT "type number of variable to be put into basis "; l PRINT IF (bas(l) <> 0) OR (min(l) < 0) OR (min(l) > -half) THEN fail2: PRINT "unreasonable variable number" GOTO fail END IF IF (l < 1) OR (l > p) THEN GOTO fail2 PRINT PRINT q = k(l) denom = a(q, l) IF denom = 0! THEN GOTO fail2 FOR i = 1 TO p IF bas(i) = q THEN bas(i) = 0 NEXT i bas(l) = q 'process tableau FOR i = 1 TO p a(q, i) = a(q, i) / denom NEXT i b(q) = b(q) / denom FOR i = 1 TO m IF i <> q THEN r = a(i, l) FOR j = 1 TO p a(i, j) = a(i, j) - r * a(q, j) NEXT j b(i) = b(i) - r * b(q) END IF NEXT i r = c(l) FOR j = 1 TO p c(j) = c(j) - r * a(q, j) NEXT j obj = obj - r * b(q) 'return for next step GOTO again END