Commit | Line | Data |
---|---|---|
7f918cf1 CE |
1 | program simplex(input, output); |
2 | ||
3 | { two-phase simplex algorithm: version Feb. 24, 1988 } | |
4 | ||
5 | { copyright K. Steiglitz } | |
6 | { Computer Science Dept. } | |
7 | { Princeton University 08544 } | |
8 | { ken@princeton.edu } | |
9 | ||
10 | const | |
11 | maxpivots = 1000; { maximum no. of pivots } | |
12 | large = 1.0e+31; { large number used in search for minimum cost column } | |
13 | lowlim = -1.0e+31; { large negative number to test for unboundedness } | |
14 | mmax = 32; { max. no. of rows } | |
15 | ncolmax = 50; { max. no. of columns allowed in tableau } | |
16 | eps = 1.0e-8; { for testing for zero } | |
17 | ||
18 | var | |
19 | done, unbounded, optimal: boolean; { flags for simplex } | |
20 | result: (toomanycols, unbound, infeas, toomanypivots, opt); | |
21 | m: 1..mmax; { no. of rows - 1, the rows are numbered 0..m } | |
22 | numpivots: integer; { pivot count } | |
23 | pivotcol, pivotrow: integer; { pivot column and row } | |
24 | pivotel: real; { pivot element } | |
25 | cbar: real; { price when searching for entering column } | |
26 | carry: array[-1..mmax, -1..mmax] of real; { inverse-basis matrix of the | |
27 | revised simplex method } | |
28 | phase: 1..2; { phase } | |
29 | price: array[0..mmax] of real; { shadow prices = row -1 of carry = | |
30 | -dual variables } | |
31 | basis: array[0..mmax] of integer; { basis columns, negative integers | |
32 | artificial } | |
33 | ncol: 1..ncolmax; { number of columns } | |
34 | tab: array[0..mmax, 1..ncolmax] of real; { tableau } | |
35 | lhs: array[0..mmax] of real; { left-hand-side } | |
36 | d: array[1..ncolmax] of real; { current cost vector } | |
37 | c: array[1..ncolmax] of real; { cost vector in original problem } | |
38 | curcol: array[-1..mmax] of real; { current column } | |
39 | curcost: real; { current cost } | |
40 | i, col, row: integer; { miscellaneous variables } | |
41 | ||
42 | procedure columnsearch; | |
43 | { looks for favorable column to enter basis. | |
44 | returns lowest cost and its column number, or turns on the flag optimal } | |
45 | ||
46 | var | |
47 | i , col : integer; | |
48 | tempcost: real; { minimum cost, temporary cost of column } | |
49 | ||
50 | begin { columnsearch } | |
51 | for i:= 0 to m do price[i]:= -carry[-1, i]; { set up price vector } | |
52 | optimal:= false; | |
53 | cbar:= large; | |
54 | pivotcol:= 0; | |
55 | for col:= 1 to ncol do | |
56 | begin | |
57 | tempcost:= d[col]; | |
58 | for i:= 0 to m do tempcost:= tempcost - price[i]*tab[i, col]; | |
59 | if( cbar > tempcost ) then | |
60 | begin | |
61 | cbar:= tempcost; | |
62 | pivotcol:= col | |
63 | end | |
64 | end; { for col } | |
65 | if ( cbar > -eps ) then optimal:= true | |
66 | end; { columnsearch } | |
67 | ||
68 | ||
69 | procedure rowsearch; | |
70 | { looks for pivot row. returns pivot row number, | |
71 | or turns on the flag unbounded } | |
72 | ||
73 | var | |
74 | i, j: integer; | |
75 | ratio, minratio: real; { ratio and minimum ratio for ratio test } | |
76 | ||
77 | begin { rowsearch } | |
78 | for i:= 0 to m do { generate column } | |
79 | begin | |
80 | curcol[i]:= 0.0; { current column = B inverse * original col. } | |
81 | for j:= 0 to m do curcol[i]:= | |
82 | curcol[i] + carry[i, j]*tab[j, pivotcol] | |
83 | end; | |
84 | curcol[-1]:= cbar; { first element in current column } | |
85 | pivotrow:= -1; | |
86 | minratio:= large; | |
87 | for i:= 0 to m do { ratio test } | |
88 | begin | |
89 | if( curcol[i] > eps ) then | |
90 | begin | |
91 | ratio:= carry[i, -1]/curcol[i]; | |
92 | if( minratio > ratio ) then { favorable row } | |
93 | begin | |
94 | minratio:= ratio; | |
95 | pivotrow:= i; | |
96 | pivotel:= curcol[i] | |
97 | end | |
98 | else { break tie with max pivot } | |
99 | if ( (minratio = ratio) and (pivotel < curcol[i]) ) then | |
100 | begin | |
101 | pivotrow:= i; | |
102 | pivotel:= curcol[i] | |
103 | end | |
104 | end { curcol > eps } | |
105 | end; { for i } | |
106 | if ( pivotrow = -1 ) then unbounded:= true { nothing found } | |
107 | else unbounded:= false | |
108 | end; { rowsearch } | |
109 | ||
110 | ||
111 | procedure pivot; | |
112 | { pivots } | |
113 | ||
114 | var | |
115 | i, j: integer; | |
116 | ||
117 | begin { pivot } | |
118 | basis[pivotrow]:= pivotcol; | |
119 | for j:= -1 to m do carry[pivotrow, j]:= carry[pivotrow, j]/pivotel; | |
120 | for i:= -1 to m do | |
121 | if( i<> pivotrow ) then | |
122 | for j:= -1 to m do | |
123 | carry[i, j]:= carry[i, j] - carry[pivotrow, j]*curcol[i]; | |
124 | curcost:= -carry[-1, -1] | |
125 | end; { pivot } | |
126 | ||
127 | ||
128 | procedure changephase; | |
129 | { changes phase from 1 to 2, by switching to original cost vector } | |
130 | ||
131 | var | |
132 | i, j, b: integer; | |
133 | ||
134 | begin { changephase } | |
135 | phase:= 2; | |
136 | for i:= 0 to m do if( basis[i] <= 0 ) then | |
137 | writeln( '...artificial basis element ', basis[i]:5, | |
138 | ' remains in basis after phase 1'); | |
139 | for j:= 1 to ncol do d[j]:= c[j]; { switch to original cost vector } | |
140 | for j:= -1 to m do | |
141 | begin | |
142 | carry[-1, j]:= 0.0; | |
143 | for i:= 0 to m do | |
144 | begin | |
145 | b:= basis[i]; { ignore artificial basis elements that are } | |
146 | if( b >= 1 ) then { still in basis } | |
147 | carry[-1, j]:= carry[-1, j] - c[b]*carry[i,j] | |
148 | end { for i } | |
149 | end; { for j } | |
150 | curcost:= -carry[-1, -1] | |
151 | end; { changephase } | |
152 | ||
153 | procedure setup; | |
154 | { sets up test problem, lhs = tab*x, x >= 0, min c*x } | |
155 | { nrow = number of rows; ncol = number of cols } | |
156 | { tab = tableau; lhs = constants } | |
157 | ||
158 | var | |
159 | i, j, nrow: integer; | |
160 | ||
161 | begin { setup } | |
162 | readln(nrow); { read number of rows } | |
163 | readln(ncol); { read number of columns } | |
164 | m:= nrow - 1; { rows are numbered 0..m } | |
165 | for j:= 1 to ncol do | |
166 | read(c[j]); { cost vector } | |
167 | for i:= 0 to m do | |
168 | begin | |
169 | read(lhs[i]); { left-hand-side } | |
170 | for j:= 1 to ncol do | |
171 | read(tab[i, j]) { tableau } | |
172 | end; | |
173 | ||
174 | done:= false; { initialize carry matrix, etc. } | |
175 | phase:= 1; | |
176 | for i:= -1 to m do for j:= -1 to mmax do carry[i, j]:= 0.0; | |
177 | for i:= 0 to m do carry[i, i]:= 1.0; { artificial basis } | |
178 | for i:= 0 to m do | |
179 | begin | |
180 | carry[i, -1]:= lhs[i]; { -1 col of carry = left-hand-side } | |
181 | carry[-1, -1]:= carry[-1, -1] - lhs[i] { - initial cost } | |
182 | end; | |
183 | curcost:= -carry[-1, -1]; | |
184 | for i:= 0 to m do basis[i]:= -i; { initial, artificial basis } | |
185 | if( ncol <= ncolmax ) then { check number of columns } | |
186 | for col:= 1 to ncol do { initialize cost vector for phase 1 } | |
187 | begin | |
188 | d[col]:= 0.0; | |
189 | for row:= 0 to m do d[col]:= d[col] - tab[row, col] | |
190 | end | |
191 | else | |
192 | begin | |
193 | writeln('...termination: too many columns for storage'); | |
194 | done:= true; | |
195 | result:= toomanycols | |
196 | end; | |
197 | numpivots:= 0; | |
198 | end; { setup } | |
199 | ||
200 | ||
201 | begin { simplex } | |
202 | setup; | |
203 | while( (numpivots < maxpivots) and (not done) and | |
204 | ( (curcost > lowlim) or (phase = 1) ) ) do | |
205 | begin | |
206 | columnsearch; | |
207 | if( not optimal ) then | |
208 | begin { not optimal } | |
209 | rowsearch; | |
210 | if( unbounded ) then | |
211 | begin | |
212 | done:= true; | |
213 | result:= unbound; | |
214 | writeln('problem is unbounded') | |
215 | end | |
216 | else | |
217 | begin | |
218 | pivot; | |
219 | numpivots:= numpivots + 1; | |
220 | if ( (numpivots = 1 ) or ( numpivots mod 10 = 0 ) ) then | |
221 | writeln('pivot ', numpivots:4, ' cost= ', curcost:12) | |
222 | end | |
223 | end { not optimal } | |
224 | else { optimal } | |
225 | if( phase = 1 ) then | |
226 | begin | |
227 | if( curcost > eps ) then | |
228 | begin | |
229 | done:= true; | |
230 | result:= infeas; | |
231 | writeln('problem is infeasible') | |
232 | end | |
233 | else | |
234 | begin | |
235 | if ( (numpivots <> 1 ) and ( numpivots mod 10 <> 0 ) ) then | |
236 | writeln('pivot ', numpivots:4, ' cost= ', curcost:12); | |
237 | writeln('phase 1 successfully completed'); | |
238 | changephase | |
239 | end | |
240 | end { if phase = 1 } | |
241 | else | |
242 | begin | |
243 | if ( (numpivots <> 1 ) and ( numpivots mod 10 <> 0 ) ) then | |
244 | writeln('pivot ', numpivots:4, ' cost= ', curcost:12); | |
245 | writeln('phase 2 successfully completed'); | |
246 | done:= true; | |
247 | result:= opt | |
248 | end | |
249 | end; { while } | |
250 | if( (curcost <= lowlim) and (phase = 2) ) then | |
251 | begin | |
252 | if ( (numpivots <> 1 ) and ( numpivots mod 10 <> 0 ) ) then | |
253 | writeln('pivot ', numpivots:4, ' cost= ', curcost:12); | |
254 | result:= unbound; | |
255 | writeln('problem is unbounded') | |
256 | end; | |
257 | if ( numpivots >= maxpivots ) then | |
258 | begin | |
259 | writeln('...termination: maximum number of pivots exceeded'); | |
260 | result:= toomanypivots | |
261 | end; | |
262 | ||
263 | if result = opt then | |
264 | begin | |
265 | writeln('optimal solution reached'); | |
266 | writeln('cost =', -carry[-1,-1]:10:6); | |
267 | for i:= 0 to m do | |
268 | writeln('x(', basis[i]:4, ')= ', carry[i,-1]:10:6) | |
269 | end | |
270 | ||
271 | end. |