'**************************************************************** '* LINEAR PROGRAMMING: THE SIMPLEX METHOD * '* ------------------------------------------------------------ * '* SAMPLE RUN: * '* Maximize z = x1 + x2 + 3x3 -0.5x4 with conditions: * '* x1 + 2x3 <= 740 * '* 2x2 - 7x4 <= 0 * '* x2 - x3 + 2x4 >= 0.5 * '* x1 + x2 + x3 +x4 = 9 * '* and all x's >=0. * '* * '* Number of variables in E.F.: 4 * '* Number of <= inequalities..: 2 * '* Number of >= inequalities..: 1 * '* Number of = equalities.....: 1 * '* Input Economic Function: * '* Coefficient # 1: 1 * '* Coefficient # 2: 1 * '* Coefficient # 3: 3 * '* Coefficient # 4: -0.5 * '* Constant term..: 0 * '* Input constraint # 1: * '* Coefficient # 1: 1 * '* Coefficient # 2: 0 * '* Coefficient # 3: 2 * '* Coefficient # 4: 0 * '* Constant term..: 740 * '* Input constraint # 2: * '* Coefficient # 1: 0 * '* Coefficient # 2: 2 * '* Coefficient # 3: 0 * '* Coefficient # 4: -7 * '* Constant term..: 0 * '* Input constraint # 3: * '* Coefficient # 1: 0 * '* Coefficient # 2: 1 * '* Coefficient # 3: -1 * '* Coefficient # 4: 2 * '* Constant term..: 0.5 * '* Input constraint # 4: * '* Coefficient # 1: 1 * '* Coefficient # 2: 1 * '* Coefficient # 3: 1 * '* Coefficient # 4: 1 * '* Constant term..: 9 * '* * '* Input Table: * '* 0.00 1.00 1.00 3.00 -0.50 * '* 740.00 -1.00 0.00 -2.00 0.00 * '* 0.00 0.00 -2.00 0.00 7.00 * '* 0.50 0.00 -1.00 1.00 -2.00 * '* 9.00 -1.00 -1.00 -1.00 -1.00 * '* * '* Maximum of E.F. = 17.02500 * '* X 1 = 0.000000 * '* X 2 = 3.325000 * '* X 3 = 4.725000 * '* X 4 = 0.950000 * '* * '* ------------------------------------------------------------ * '* Reference: "Numerical Recipes By W.H. Press, B. P. Flannery, * '* S.A. Teukolsky and W.T. Vetterling, Cambridge * '* University Press, 1986" [BIBLI 08]. * '* * '* Basic Release 1.0 By J-P Moreau, Paris * '* (www.jpmoreau.fr) * '**************************************************************** 'Program Test Simplex DEFINT I-N MMAX = 25: NMAX = 25 DIM A(MMAX, NMAX) DIM iposv(MMAX), izrov(NMAX) CLS PRINT PRINT " Number of variables in E.F.: "; : INPUT "", N PRINT " Number of <= inequalities..: "; : INPUT "", M1 PRINT " Number of >= inequalities..: "; : INPUT "", M2 PRINT " Number of = equalities.....: "; : INPUT "", M3 M = M1 + M2 + M3 'Total number of constraints F$ = "#####.##" F1$ = "#####.######" FOR i = 1 TO M + 2 FOR j = 1 TO N + 1 A(i, j) = 0! NEXT j NEXT i PRINT " Input Economic Function:" FOR i = 2 TO N + 1 PRINT " Coefficient #"; i - 1; ": "; INPUT "", A(1, i) NEXT i PRINT " Constant term : "; INPUT "", A(1, 1) ' input constraints FOR i = 1 TO M PRINT " Input constraint #"; i; ":" FOR j = 2 TO N + 1 PRINT " Coefficient #"; j - 1; ": "; INPUT "", R A(i + 1, j) = -R NEXT j PRINT " Constant term : "; INPUT "", A(i + 1, 1) NEXT i PRINT PRINT " Input Table:" FOR i = 1 TO M + 1 FOR j = 1 TO N + 1 PRINT USING F$; A(i, j); NEXT j PRINT NEXT i GOSUB 1000 'call simplx(A,M,N,M1,M2,M3,ICASE,IZROV,IPOSV) IF icase = 0 THEN 'result ok. PRINT PRINT " Maximum of E.F. = "; PRINT USING F1$; A(1, 1) FOR i = 1 TO N FOR j = 1 TO M IF iposv(j) = i THEN PRINT " X"; i; " = "; PRINT USING F1$; A(j + 1, 1) GOTO 3 END IF NEXT j PRINT " X"; i; " = "; PRINT USING F1$; 0! 3 NEXT i ELSE PRINT " No solution (error code = "; icase; ")." END IF PRINT END 'of main program 1000 'Subroutine simplx(A, m, n, m1, m2, m3, icase, izrov, iposv) '----------------------------------------------------------------------------------------- 'USES 2000 simp1, 2100 simp2, 2200 simp3. 'Simplex method for linear programming. Input parameters A, m, n, m1, m2, and m3, and 'output parameters A, icase, izrov, and iposv are described above (see reference). 'Constants: MMAX is the maximum number of constraints expected; NMAX is the maximum number 'of variables expected; EPS is the absolute precision, which should be adjusted to the 'scale of your variables. '----------------------------------------------------------------------------------------- 'Labels: 1,2,10,20,30 DIM l1(MMAX), l2(MMAX), l3(MMAX) EPS = .000001 IF M <> M1 + M2 + M3 THEN PRINT " Bad input constraint counts in simplx." RETURN END IF nl1 = N FOR k = 1 TO N l1(k) = k 'Initialize index list of columns admissible for exchange. izrov(k) = k 'Initially make all variables right-hand. NEXT k nl2 = M FOR i = 1 TO M IF A(i + 1, 1) < 0! THEN PRINT " Bad input tableau in simplx, Constants bi must be nonnegative." RETURN END IF l2(i) = i iposv(i) = N + i '------------------------------------------------------------------------------------------------- 'Initial left-hand variables. m1 type constraints are represented by having their slack variable 'initially left-hand, with no artificial variable. m2 type constraints have their slack 'variable initially left-hand, with a minus sign, and their artificial variable handled implicitly 'during their first exchange. m3 type constraints have their artificial variable initially 'left-hand. '------------------------------------------------------------------------------------------------- NEXT i FOR i = 1 TO M2 l3(i) = 1 NEXT i ir = 0 IF M2 + M3 = 0 THEN GOTO 30 'The origin is a feasible starting solution. Go to phase two. ir = 1 FOR k = 1 TO N + 1 'Compute the auxiliary objective function. q1 = 0! FOR i = M1 + 1 TO M q1 = q1 + A(i + 1, k) NEXT i A(M + 2, k) = -q1 NEXT k 10 : iabf = 0: mm = M + 1 GOSUB 2000 'call simp1(a,m+1,l1,nl1,0,kp,bmax) Find max. coeff. of auxiliary objective fn. IF bmax <= EPS AND A(M + 2, 1) < -EPS THEN icase = -1 'Auxiliary objective function is still negative and can’t be improved, RETURN 'hence no feasible solution exists. ELSEIF bmax <= EPS AND A(M + 2, 1) <= EPS THEN ' Auxiliary objective function is zero and can’t be improved; we have a feasible starting vector. ' Clean out the artificial variables corresponding to any remaining equality constraints by ' go to 1’s and then move on to phase two by go to 30. m12 = M1 + M2 + 1 IF m12 <= M THEN FOR ip = m12 TO M IF iposv(ip) = ip + N THEN 'Found an artificial variable for an equality constraint. iabf = 1: mm = ip GOSUB 2000 'call simp1(a,ip,l1,nl1,1,kp,bmax) IF bmax > EPS THEN GOTO 1 'Exchange with column corresponding to maximum END IF 'pivot element in row. NEXT ip END IF ir = 0 m12 = m12 - 1 IF M1 + 1 > m12 THEN GOTO 30 FOR i = M1 + 1 TO M1 + M2 'Change sign of row for any m2 constraints IF l3(i - M1) = 1 THEN 'still present from the initial basis. FOR k = 1 TO N + 1 A(i + 1, k) = -1! * A(i + 1, k) NEXT k END IF NEXT i GOTO 30 'Go to phase two. END IF GOSUB 2100 'call simp2(a,m,n,l2,nl2,ip,kp,q1) Locate a pivot element (phase one). IF ip = 0 THEN 'Maximum of auxiliary objective function is icase = -1 'unbounded, so no feasible solution exists. RETURN END IF 1 i1 = M + 1: k1 = N GOSUB 2200 'call simp3(a,m+1,n,ip,kp) ' Exchange a left- and a right-hand variable (phase one), then update lists. IF iposv(ip) >= N + M1 + M2 + 1 THEN 'Exchanged out an artificial variable for an 'equality constraint. Make sure it stays FOR k = 1 TO nl1 'out by removing it from the l1 list. IF l1(k) = kp THEN GOTO 2 NEXT k 2 nl1 = nl1 - 1 FOR iis = k TO nl1 l1(iis) = l1(iis + 1) NEXT iis ELSE IF iposv(ip) < N + M1 + 1 THEN GOTO 20 kh = iposv(ip) - M1 - N IF l3(kh) = 0! THEN GOTO 20 'Exchanged out an m2 type constraint. l3(kh) = 0! 'If it’s the first time, correct the pivot column or the 'minus sign and the implicit artificial variable. END IF A(M + 2, kp + 1) = A(M + 2, kp + 1) + 1! FOR i = 1 TO M + 2 A(i, kp + 1) = -1! * A(i, kp + 1) NEXT i 20 iis = izrov(kp) 'Update lists of left- and right-hand variables. izrov(kp) = iposv(ip) iposv(ip) = iis IF ir <> 0 THEN GOTO 10 'if still in phase one, go back to 10. ' End of phase one code for finding an initial feasible solution. Now, in phase two, optimize it. 30 : iabf = 0: mm = 0 GOSUB 2000 'call simp1(a,0,l1,nl1,0,kp,bmax) Test the z-row for doneness. IF bmax <= EPS THEN 'Done. Solution found. Return with the good news. icase = 0 RETURN END IF GOSUB 2100 'call simp2(a,m,n,l2,nl2,ip,kp,q1) Locate a pivot element (phase two). IF ip = 0 THEN 'Objective function is unbounded. Report and return. icase = 1 RETURN END IF i1 = M: k1 = N GOSUB 2200 'call simp3(a,m,n,ip,kp) Exchange a left- and a right-hand variable (phase two), GOTO 20 'update lists of left- and right-hand variables and 'return for another iteration. RETURN 'The preceding subroutine makes use of the following utility subroutines: 2000 'Subroutine simp1(a, mm, l1, nl1, iabf, kp, bmax) ' Determines the maximum of those elements whose index is contained in the supplied list ' l1, either with or without taking the absolute value, as flagged by iabf. kp = l1(1) bmax = A(mm + 1, kp + 1) IF nl1 < 2 THEN RETURN FOR k = 2 TO nl1 IF iabf = 0 THEN test = A(mm + 1, l1(k) + 1) - bmax ELSE test = ABS(A(mm + 1, l1(k) + 1)) - ABS(bmax) END IF IF test > 0! THEN bmax = A(mm + 1, l1(k) + 1) kp = l1(k) END IF NEXT k RETURN 2100 'Subroutine simp2(a, m, n, l2, nl2, ip, kp, q1) EPS = .000001 ' Locate a pivot element, taking degeneracy into account. ip = 0 IF nl2 < 1 THEN RETURN FOR i = 1 TO nl2 IF A(i + 1, kp + 1) < -EPS THEN GOTO 2102 NEXT i RETURN 'No possible pivots. Return with message. 2102 q1 = -A(l2(i) + 1, 1) / A(l2(i) + 1, kp + 1) ip = l2(i) IF i + 1 > nl2 THEN RETURN FOR i = i + 1 TO nl2 ii = l2(i) IF A(ii + 1, kp + 1) < -EPS THEN q = -A(ii + 1, 1) / A(ii + 1, kp + 1) IF q < q1 THEN ip = ii q1 = q ELSEIF q = q1 THEN 'We have a degeneracy. FOR k = 1 TO N qp = -A(ip + 1, k + 1) / A(ip + 1, kp + 1) q0 = -A(ii + 1, k + 1) / A(ii + 1, kp + 1) IF q0 <> qp THEN GOTO 2106 NEXT k 2106 IF q0 < qp THEN ip = ii END IF END IF NEXT i RETURN 2200 'Subroutine simp3(a, i1,k1,ip,kp) ' Matrix operations to exchange a left-hand and right-hand variable. piv = 1! / A(ip + 1, kp + 1) IF i1 >= 0 THEN FOR ii = 1 TO i1 + 1 IF ii - 1 <> ip THEN A(ii, kp + 1) = A(ii, kp + 1) * piv FOR kk = 1 TO k1 + 1 IF kk - 1 <> kp THEN A(ii, kk) = A(ii, kk) - A(ip + 1, kk) * A(ii, kp + 1) END IF NEXT kk END IF NEXT ii END IF FOR kk = 1 TO k1 + 1 IF kk - 1 <> kp THEN A(ip + 1, kk) = -A(ip + 1, kk) * piv NEXT kk A(ip + 1, kp + 1) = piv RETURN ' end of file tsimplex.bas