← Back to branch summary

~ubuntu-branches/debian/squeeze/maxima/squeeze

« back to all changes in this revision

Viewing changes to share/contrib/simplex/simplex_algorithm.lisp

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2006-10-18 14:52:42 UTC
  • mto: (1.1.5 upstream)
  • mto: This revision was merged to the branch mainline in revision 4.
  • Revision ID: james.westby@ubuntu.com-20061018145242-vzyrm5hmxr8kiosf
Import upstream version 5.10.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
share/contrib/simplex/simplex_algorithm.lisp
 
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_sx:     epsilon for numeric computation (default: 1e-8)         ;;
 
79
;; - scale_sx:       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
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
85
 
 
86
(in-package :maxima)
 
87
(macsyma-module simplex)
 
88
 
 
89
(proclaim '(optimize (speed 3) (safety 0) (space 0) (compilation-speed 0)))
 
90
 
 
91
($put '$simplex 1.02 '$version)
 
92
 
 
93
(defmvar $pivot_count_sx     0  "Number of pivots in last problem."   fixnum)
 
94
(defmvar $pivot_max_sx   15000  "Maximum number of pivots allowed."   fixnum)
 
95
(defmvar $epsilon_sx      1e-8  "Epsilon for numerical computation."  flonum)
 
96
(defmvar $scale_sx         nil  "Should we scale the input."         boolean)
 
97
 
 
98
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
99
;;                                                                           ;;
 
100
;; Two-phase standard simplex method for solving linear program in standard  ;;
 
101
;; form.                                                                     ;;
 
102
;;                                                                           ;;
 
103
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
104
 
 
105
(defun $linear_program (A b c)
 
106
  (if (not ($matrixp A))
 
107
      (merror "linear_program: first argument not matrix."))
 
108
  (if (not ($listp b))
 
109
      (merror "linear_program: second argument not list."))
 
110
  (if (not ($listp c))
 
111
      (merror "linear_program: third argument not list."))
 
112
  (if (not (meqp ($length b) ($length A)))
 
113
      (merror "linear_program: second argument not of correct length."))
 
114
  (if (not (meqp ($length c) ($length ($first A))))
 
115
      (merror "linear_program: third argument not of correct length."))
 
116
  
 
117
  (let* ((m ($length A))
 
118
         (n ($length ($first A)))
 
119
         (Tab (make-array `(,(+ 2 m) ,(1+ n)))) ; Tableau
 
120
         (basis ())      ; which columns are in current basis
 
121
         (sc-fac ()))    ; scaling factors
 
122
    
 
123
    (setq $pivot_count_sx 0)
 
124
 
 
125
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
126
;; Construct the tableau for phase 1:                              ;;
 
127
;;       [   A      b    ]                                         ;;
 
128
;; Tab = [   c'     0    ]                                         ;;
 
129
;;       [ sm(A)  sm(b)  ]                                         ;;
 
130
;;                                                                 ;;
 
131
;; If b[i]<0 multiply b[i] and A[i] by -1 (required for phase 1).  ;;
 
132
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
133
    (dotimes (i m)
 
134
      (if (mlsp (maref b (1+ i)) 0)
 
135
          (progn
 
136
            (dotimes (j n)
 
137
              (setf (aref Tab i j) (neg (maref A (1+ i) (1+ j)))))
 
138
            (setf (aref Tab i n) (neg (maref b (1+ i)))))
 
139
          (progn
 
140
            (dotimes (j n)
 
141
              (setf (aref Tab i j) (maref A (1+ i) (1+ j))))
 
142
            (setf (aref Tab i n) (maref b (1+ i))))))
 
143
    (dotimes (i n)
 
144
      (setf (aref Tab m i) (neg (maref c (1+ i)))))
 
145
    (dotimes (i n)
 
146
      (setf (aref Tab (1+ m) i) 0)
 
147
      (dotimes (j m)
 
148
        (setf (aref Tab (1+ m) i) (add (aref Tab (1+ m) i)
 
149
                                       (aref Tab j i)))))
 
150
    (setf (aref Tab (1+ m) n) 0)
 
151
    (dotimes (i m)
 
152
      (setf (aref Tab (1+ m) n) (add (aref Tab (1+ m) n)
 
153
                                     (aref Tab i n))))
 
154
    (setf (aref Tab m n) 0)
 
155
 
 
156
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
157
;; At the beginning the artificial variables are in the basis.     ;;
 
158
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
159
    (dotimes (i m)
 
160
      (setq basis (append basis (list (add n i)))))
 
161
 
 
162
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
163
;; Scaling.                                                        ;;
 
164
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
165
    (dotimes (i n)
 
166
      (setq sc-fac (append sc-fac (list 1))))
 
167
    (if $scale_sx
 
168
        (scale-sx Tab (+ 2 m) (1+ n) sc-fac))
 
169
 
 
170
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
171
;; Phase 1                                                         ;;
 
172
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
173
 
 
174
    (simplex-sx Tab basis m (+ 2 m) (1+ n))
 
175
 
 
176
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
177
 
 
178
    (if (mlsp $epsilon_sx (aref Tab (1+ m) n))
 
179
        "Problem not feasible!"
 
180
        (let ((is-bounded))
 
181
 
 
182
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
183
;; Check for artificial variables in basis.                        ;;
 
184
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
185
 
 
186
          (dotimes (i m)
 
187
            (if (>= (nth i basis) n)
 
188
                (if (not (run-out-of-basis-sx Tab m n basis i))
 
189
                    ($print "Matrix A is not of full rank:"
 
190
                            "Row" (1+ i) "is redundant!"))))
 
191
 
 
192
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
193
;; Phase 2                                                         ;;
 
194
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
195
 
 
196
          (setq is-bounded (simplex-sx Tab basis m (1+ m) (1+ n)))
 
197
 
 
198
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
199
 
 
200
          (if is-bounded
 
201
              (let ((opt-tmp ())
 
202
                    (opt ()))
 
203
                (dotimes (i (+ m n))
 
204
                  (setq opt-tmp (append opt-tmp (list 0))))
 
205
                (dotimes (i m)
 
206
                  (setf (nth (nth i basis) opt-tmp)
 
207
                        (div (aref Tab i n)                ;; undo the
 
208
                             (nth (nth i basis) sc-fac)))) ;; scaling
 
209
                (dotimes (i n)
 
210
                  (setq opt (append opt (list 0))))
 
211
                (dotimes (i n)
 
212
                  (setf (nth i opt) (nth i opt-tmp)))
 
213
                (setq opt (cons '(mlist simp) opt))        ;; build the
 
214
                (setq opt `((mlist simp) ,opt              ;; the solution
 
215
                            ,(aref Tab m n))))
 
216
              "Problem not bounded!")))))
 
217
 
 
218
 
 
219
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
220
;;                                                                           ;;
 
221
;; Simplex algorithm.                                                        ;;
 
222
;;                                                                           ;;
 
223
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
224
 
 
225
(defun simplex-sx (Tab basis Am m n)
 
226
  (let ((ip) (jp) (tmp)  (is-bounded))
 
227
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
228
  ;; Repeat while we don't have a solution or know that the        ;;
 
229
  ;; problem is unbounded.                                         ;;
 
230
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
231
    (do ((have-solution nil))
 
232
        ((not (null have-solution)))
 
233
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
234
  ;; Choose jp so that jp-th column has the biggest reduced cost.  ;;
 
235
  ;; If all reduced costs are negative, we have the solution.      ;;
 
236
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
237
      (setq tmp (aref Tab (1- m) 0))
 
238
      (setq jp 0)
 
239
      (dotimes (j (- n 1))
 
240
        (if (mlsp tmp (aref Tab (1- m) j))
 
241
            (progn
 
242
              (setq tmp (aref Tab (1- m) j))
 
243
              (setq jp j))))
 
244
      (if (or (mlsp tmp $epsilon_sx)
 
245
              (meqp tmp $epsilon_sx))
 
246
          (progn
 
247
            (setq is-bounded t)
 
248
            (setq have-solution t))
 
249
          (progn
 
250
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
251
  ;; Choose ip so that Tab[ip,n]/Tab[ip,jp] is the smallest        ;;
 
252
  ;; possible among those for which Tab[i,n] is positive. If all   ;;
 
253
  ;; Tab[i,n] are negative, the problem is unbounded.              ;;
 
254
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
255
            (setq tmp nil)
 
256
            (setq ip 0)
 
257
            (dotimes (i Am)
 
258
              (if (mlsp $epsilon_sx (aref Tab i jp))
 
259
                  (if (or (null tmp) (mlsp (div (aref Tab i (1- n))
 
260
                                                (aref Tab i jp))
 
261
                                           tmp))
 
262
                      (progn
 
263
                        (setq tmp (div (aref Tab i (1- n))
 
264
                                       (aref Tab i jp)))
 
265
                        (setq ip i)))))
 
266
            (if (null tmp)
 
267
                (progn
 
268
                  (setq is-bounded nil)
 
269
                  (setq have-solution t))
 
270
                (progn
 
271
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
272
  ;; Pivot the simplex tableau.                                   ;;
 
273
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
274
                  (setf (nth ip basis) jp)
 
275
                  (pivot-sx Tab ip jp m n))))))
 
276
    is-bounded))
 
277
 
 
278
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
279
;;                                                                           ;;
 
280
;; Pivoting for the Simplex algorithm.                                       ;;
 
281
;;                                                                           ;;
 
282
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
283
 
 
284
(defun pivot-sx (Tab ip jp m n)
 
285
  (let ((piv (aref Tab ip jp)))
 
286
    (setq $pivot_count_sx (1+ $pivot_count_sx))
 
287
    (if (meqp $pivot_count_sx $pivot_max_sx)
 
288
        (progn
 
289
          ($print "Maximum number of pivots reached.")
 
290
          ($print "Try setting a bigger value for pivot_max_sx.")
 
291
          ($print "Try setting scale_sx to true.")
 
292
          ($print "")
 
293
          ($error "linear_program: maximum number of pivots reached.")))
 
294
    (if (meqp piv 0)
 
295
        ($print "Singular!")
 
296
        (progn
 
297
          (dotimes (i n)
 
298
            (setf (aref Tab ip i) (div (aref Tab ip i) piv)))
 
299
          (dotimes (i m)
 
300
            (if (not (eq i ip))
 
301
                (let ((tm (aref Tab i jp)))
 
302
                  (if (not (meqp tm 0))
 
303
                      (dotimes (j n)
 
304
                        (setf (aref Tab i j)
 
305
                              (sub (aref Tab i j)
 
306
                                   (mul tm (aref Tab ip j)))))))))))))
 
307
 
 
308
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
309
;;                                                                           ;;
 
310
;; Run artificial variable out of basis.                                     ;;
 
311
;;                                                                           ;;
 
312
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
313
 
 
314
(defun run-out-of-basis-sx (Tab m n basis i)
 
315
  (let ((jp nil))
 
316
    (do ((j 0 (1+ j)))
 
317
        ((or (= j n) (not (null jp))))
 
318
      (if (not (meqp 0 (aref Tab i j))) ; if Tab[i,j]#0 then column j is not
 
319
          (setq jp j)))                 ; in the basis
 
320
    (if (null jp)
 
321
        nil
 
322
        (progn
 
323
          (setf (nth i basis) jp)
 
324
          (pivot-sx Tab i jp (+ 2 m) (+ 1 n))
 
325
          1))))
 
326
 
 
327
 
 
328
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
329
;;                                                                           ;;
 
330
;; Scaling for the simplex algorithm. (Equilibratium scaling)                ;;
 
331
;;                                                                           ;;
 
332
;; After scaling, the maximum absolute value in each row/column is 1.        ;;
 
333
;;                                                                           ;;
 
334
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
335
 
 
336
(defun scale-sx (Tab m n sc-fac)
 
337
  (let ((r 0))
 
338
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
339
  ;; Scale the rows of A and b.                                   ;;
 
340
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
341
    (dotimes (i (- m 2))
 
342
      (setq r 0)
 
343
      (dotimes (j n)
 
344
        (let* ((tij (aref Tab i j))
 
345
               (ta (if (mlsp tij 0) (neg tij) tij)))
 
346
          (if (mlsp r ta)
 
347
              (setq r ta))))
 
348
      (if (mlsp $epsilon_sx r)
 
349
          (dotimes (j n)
 
350
            (setf (aref Tab i j) (div (aref Tab i j) r)))))
 
351
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
352
  ;; Scale the columns of A and c.                                ;;
 
353
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
354
    (dotimes (j (1- n))
 
355
      (setq r 0)
 
356
      (dotimes (i (- m 2))
 
357
        (let* ((tij (aref Tab i j))
 
358
               (ta (if (mlsp tij 0) (neg tij) tij)))
 
359
          (if (mlsp r ta)
 
360
              (setq r ta))))
 
361
      (if (mlsp $epsilon_sx r)
 
362
          (progn
 
363
            (dotimes (i m)
 
364
              (setf (aref Tab i j) (div (aref Tab i j) r)))
 
365
            (setf (nth j sc-fac) r))))))