| blob | 832c31f7ce56174316d6615e8fa2b6db52776a80 |
1 ;;; -*- Mode: Lisp -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;; ;;
4 ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;
5 ;; ;;; ;;; ;;
6 ;; ;;; ~*~ SIMPLEX ~*~ ;;; ;;
7 ;; ;;; ;;; ;;
8 ;; ;;; A simple implementation of the simplex ;;; ;;
9 ;; ;;; algorithm for Linear Programming for Maxima. ;;; ;;
10 ;; ;;; ;;; ;;
11 ;; ;;; ;;; ;;
12 ;; ;;; Copyright: Andrej Vodopivec <andrejv@users.sourceforge.net> ;;; ;;
13 ;; ;;; Version: 1.02 ;;; ;;
14 ;; ;;; License: GPL ;;; ;;
15 ;; ;;; ;;; ;;
16 ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;
17 ;; ;;
18 ;; ;;
19 ;; USAGE: ;;
20 ;; ======= ;;
21 ;; ;;
22 ;; To get optimal speed this file should be compiled. ;;
23 ;; ;;
24 ;; linear_program(A, b, c): ;;
25 ;; ;;
26 ;; - problem: - find x which minimizes c'.x with constraints A.x=b and x>=0 ;;
27 ;; ;;
28 ;; - input: - matrix A of size mxn, ;;
29 ;; - vector (list) b of length m, ;;
30 ;; - vector (list) c of length n ;;
31 ;; ;;
32 ;; - output: - [x, val] if the problem is solvable; ;;
33 ;; x is the optimal vector, val=c'.x ;;
34 ;; - "Problem not bounded" if the problem is not bounded ;;
35 ;; - "Problem not feasible" if the problem is not feasible ;;
36 ;; ;;
37 ;; - algorithm: a simple implementation of the two-phase simplex algorithm ;;
38 ;; ;;
39 ;; ;;
40 ;; DEMO: ;;
41 ;; ====== ;;
42 ;; ;;
43 ;; We would like to minimize x-2*y subject to: ;;
44 ;; ;;
45 ;; x + y >= 1 ;;
46 ;; 2*x - 3*y >= 1 ;;
47 ;; 4*x - 5*y = 6 ;;
48 ;; x, y >= 0 ;;
49 ;; ;;
50 ;; We have to introduce two slack variables for inequalities to get a ;;
51 ;; linear program in the standard form: ;;
52 ;; ;;
53 ;; x + y - s1 = 1 ;;
54 ;; 2*x - 3*y - s2 = 1 ;;
55 ;; 4*x - 5*y = 6 ;;
56 ;; x, y, s1, s2 >= 0 ;;
57 ;; ;;
58 ;; Construct parameters: ;;
59 ;; ;;
60 ;; A : matrix([1,1,-1,0],[2,-3,0,-1], [4,-5,0,0]); ;;
61 ;; b : [1,1,6]; ;;
62 ;; c : [1,-2,0,0]; ;;
63 ;; ;;
64 ;; Solution: ;;
65 ;; ;;
66 ;; linear_program(A, b, c); ;;
67 ;; => [[13/2, 4, 19/2, 0], -3/2] ;;
68 ;; ;;
69 ;; The solution is: x=13/2 and y=4 (s1=19/2, s2=0), and the value of the ;;
70 ;; minimum is x-2*y=-3/2. ;;
71 ;; ;;
72 ;; ;;
73 ;; VARIABLES: ;;
74 ;; =========== ;;
75 ;; ;;
76 ;; - pivot_count_sx: the number of pivots in last computation ;;
77 ;; - pivot_max_sx: the maximum number of pivots in computation ;;
78 ;; - epsilon_lp: epsilon for numeric computation (default: 1e-8) ;;
79 ;; - scale_lp: should maxima scale the problem: can be used in ;;
80 ;; Klee-Minty problem to speed-up computation or in some ;;
81 ;; cases to improve numerical stability (default: false); ;;
82 ;; uses equilibratium scaling ;;
83 ;; ;;
84 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
97 ;; type checkers
100 ;; comparison
110 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
111 ;; ;;
112 ;; Two-phase standard simplex method for solving linear program in standard ;;
113 ;; form. ;;
114 ;; ;;
115 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
139 (mwarning (format nil "linear_program(A,b,c): non-rat inputs found, epsilon_lp=~e." $epsilon_lp))
144 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
145 ;; Construct the tableau for phase 1: ;;
146 ;; [ A b ] ;;
147 ;; Tab = [ c' 0 ] ;;
148 ;; [ sm(A) sm(b) ] ;;
149 ;; ;;
150 ;; If b[i]<0 multiply b[i] and A[i] by -1 (required for phase 1). ;;
151 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
175 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
176 ;; At the beginning the artificial variables are in the basis. ;;
177 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
181 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
182 ;; Scaling. ;;
183 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
189 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
190 ;; Phase 1 ;;
191 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
195 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
198 "Problem not feasible!"
201 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
202 ;; Check for artificial variables in basis. ;;
203 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
208 $warn_rank_sx)
212 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
213 ;; Phase 2 ;;
214 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
218 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
239 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
240 ;; ;;
241 ;; Simplex algorithm. ;;
242 ;; ;;
243 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
247 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
248 ;; Repeat while we don't have a solution or know that the ;;
249 ;; problem is unbounded. ;;
250 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
253 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
254 ;; Choose jp so that jp-th column has the biggest reduced cost. ;;
255 ;; If all reduced costs are negative, we have the solution. ;;
256 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
269 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
270 ;; Choose ip so that Tab[ip,n]/Tab[ip,jp] is the smallest ;;
271 ;; possible among those for which Tab[i,n] is positive. If all ;;
272 ;; Tab[i,n] are negative, the problem is unbounded. ;;
273 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
280 tmp))
290 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
291 ;; Pivot the simplex tableau. ;;
292 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
295 is-bounded))
297 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
298 ;; ;;
299 ;; Pivoting for the Simplex algorithm. ;;
300 ;; ;;
301 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
327 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
328 ;; ;;
329 ;; Run artificial variable out of basis. ;;
330 ;; ;;
331 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
340 nil
347 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
348 ;; ;;
349 ;; Scaling for the simplex algorithm. (Equilibratium scaling) ;;
350 ;; ;;
351 ;; After scaling, the maximum absolute value in each row/column is 1. ;;
352 ;; ;;
353 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
357 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
358 ;; Scale the rows of A and b. ;;
359 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
370 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
371 ;; Scale the columns of A and c. ;;
372 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;