home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / language / examples / simplex.lsp < prev    next >
Lisp/Scheme  |  1993-10-23  |  27KB  |  578 lines
  1. ;; Simplex-Algorithmus und dessen Demonstration
  2. ;; Bruno Haible 25.09.1986, 30.5.1991
  3. ;; Version für Common Lisp
  4.  
  5. (provide 'simplex)
  6.  
  7. ; Typ einer Zahl: hier rationale Zahlen
  8. (deftype Zahl () 'rational) ; abgekürzt: R
  9. (deftype Zahl-Vektor () '(simple-array Zahl 1)) ; abgekürzt: R^n
  10. (deftype Zahl-Matrix () '(simple-array Zahl 2)) ; abgekürzt: R^mxn
  11.  
  12. ; Input: A, eine mxn-Matrix von Zahlen,
  13. ;        b, ein m-Vektor von Zahlen,
  14. ;        c, ein n-Vektor von Zahlen.
  15. ; Löst eine kanonische Optimierungsaufgabe :
  16. ; ( C ) Finde unter den y ∈ R^n mit A*y=b, y≥0 diejenigen mit <c,y> = Min
  17. ; ( C* ) Finde unter den z ∈ R^m mit x=A`*z+c≥0 diejenigen mit <b,z> = Min
  18. ; Siehe Koulikov, Algèbre et théorie des nombres, Chap. 6.1.
  19. ; Output: 1. Wert: ob lösbar, NIL oder T,
  20. ;         falls lösbar:
  21. ;           2. Wert: Lösung von ( C ), dem LP: Liste, bestehend aus
  22. ;                    - der Lösung y
  23. ;                    - dem optimalen Wert von <c,y>
  24. ;           3. Wert: Lösung von ( C* ), dem Dualproblem DP: Liste, bestehend aus
  25. ;                    - der Lösung x
  26. ;                    - den konstanten Gliedern Z für die z's
  27. ;                    - den Abhängigkeiten ZZ zwischen den z's
  28. ;                    - dem optimalen Wert von <b,z>
  29. ;                    Die allgemeine Lösung ist so aufgebaut:
  30. ;                    Falls ZZ[i,i]=T, ist die Variable z_i frei wählbar.
  31. ;                    Sonst ist ZZ[i,i]=NIL, und die Variable z_i errechnet sich als
  32. ;                    z_i = Z[i] + Summe(j=1,...,m mit Z[j,j]=T , ZZ[i,j]*z_j) .
  33. ;           4. Wert: ob die angegebenen Lösungen die gesamte Lösungsmenge
  34. ;                    ausmachen (T) oder ob möglicherweise noch eine andere
  35. ;                    Lösung existiert.
  36. ;                    In jedem Fall ist die Lösungsmenge konvex.
  37. ; Methode:
  38. ; Methode von Dantzig, siehe Koulikov-Buch;
  39. ; Zusatz für entarteten Fall nach Bronstein/Semendjajew.
  40. (defun simplex (A b c)
  41.   ; Input überprüfen:
  42.   (let ((m (length b))
  43.         (n (length c)))
  44.     (declare (type fixnum m n))
  45.     (assert (typep A `(array * (,m ,n))))
  46.     ; m = die Höhe der Matrix, in der gearbeitet wird.
  47.     ; n = die Breite der Matrix, in der gearbeitet wird.
  48.     ; Initialisierung:
  49.     (let ((AA (let ((h (make-array `(,m ,n))))
  50.                 (dotimes (i m)
  51.                   (dotimes (j n)
  52.                     (setf (aref h i j) (rational (aref A i j)))
  53.                 ) )
  54.                 h
  55.           )   )
  56.           (AB (let ((h (make-array m)))
  57.                 (dotimes (i m) (setf (aref h i) (rational (elt b i))))
  58.                 h
  59.           )   )
  60.           (AC (let ((h (make-array n)))
  61.                 (dotimes (j n) (setf (aref h j) (rational (elt c j))))
  62.                 h
  63.           )   )
  64.           (AD 0)
  65.           ; AA ∈ R^mxn : das Tableau, in dem gearbeitet wird
  66.           ; AB ∈ R^m : der rechte Rand
  67.           ; AC ∈ R^n : der untere Rand
  68.           ; AD ∈ R : die rechte untere Ecke
  69.           (Extend nil)
  70.           (AAE nil)
  71.           (ACE nil)
  72.           ; Extend ∈ Boolean : sind die Zusatzspalten (Entartungsfall) aktiv?
  73.           ; AAE ∈ R^mxm : Zusatzspalten für den Entartungsfall
  74.           ; ACE ∈ R^m : unterer Rand der Zusatzspalten
  75.           (Zeile (let ((h (make-array m)))
  76.                    (dotimes (i m) (setf (aref h i) (cons 'ZO i)))
  77.                    h
  78.           )      )
  79.           (Spalte (let ((h (make-array n)))
  80.                     (dotimes (j n) (setf (aref h j) (cons 'XY j)))
  81.                     h
  82.           )       )
  83.           ; Zeile ∈ Cons^m : linke Markierungen
  84.           ; Spalte ∈ Cons^n : obere Markierungen
  85.           ;
  86.           ; Skizze des Tableaus:
  87.           ;
  88.           ;             | Zeile[0]  Zeile[1]  ... Zeile[n-1]  |
  89.           ; ------------+-------------------------------------+---------
  90.           ; Spalte[0]   |  AA[0,0]   AA[0,1]  ... AA[0,n-1]   |  AB[0]
  91.           ; Spalte[1]   |  AA[1,0]   AA[1,1]  ... AA[1,n-1]   |  AB[1]
  92.           ;   ...       |    ...       ...          ...       |   ...
  93.           ; Spalte[m-1] | AA[m-1,0] AA[m-1,1] ... AA[m-1,n-1] | AB[m-1]
  94.           ; ------------+-------------------------------------+---------
  95.           ;             |   AC[0]     AC[1]   ...   AC[n-1]   |   AD
  96.           ;
  97.           ; Markierung von Spalten:   0   oder   Y
  98.           ;                          ...        ...
  99.           ;                           Z          X
  100.           ;
  101.           ; Markierung von Zeilen:  Z ... 0  oder  X ... -Y.
  102.           ;
  103.           ; für alle i=0,...,m-1:
  104.           ;   Summe(j=0,...,n-1; AA[i,j]*(0 oder Y[Spalte[j].Nr])) - AB[i] =
  105.           ;     = (0 oder -Y[Zeile[i].Nr])
  106.           ; für alle j=0,...,n-1:
  107.           ;   Summe(i=0,...,m-1; AA[i,j]*(Z[Zeile[i].Nr] oder X[Zeile[i].Nr])) + AC[j] =
  108.           ;     = (Z[Spalte[j].Nr] oder X[Spalte[j].Nr])
  109.           ; Summe(j=1,...,N; AC[j]*(0 oder Y[Spalte[j].Nr])) - AD = <c,y>
  110.           ; Summe(i=1,...,M; AB[i]*(Z[Zeile[i].Nr] oder X[Zeile[i].Nr])) + AD = <b,z>
  111.           ; Dieses ist als Gleichung mit den Unbekannten X,Y,Z zu verstehen.
  112.           ;
  113.           ; Die Zusatzspalten sind - wenn vorhanden - rechts anzufügen.
  114.           ;
  115.           (ZZ (let ((h (make-array `(,m ,m) :initial-element 0)))
  116.                 (dotimes (i m) (setf (aref h i i) nil))
  117.                 h
  118.           )   )
  119.           (Z (make-array m))
  120.           (Y (make-array n))
  121.           (X (make-array n))
  122.          )
  123.       (declare (type Zahl-Matrix AA) (type Zahl-Vektor AB AC) (type Zahl AD))
  124.       (flet
  125.         ((pivot (k l)
  126.            (declare (type fixnum k l))
  127.            ; pivotiert das Tableau das Element AA[k,l] mit 0≤k<m, 0≤l<n.
  128.            ; Invariante ist, daß vor und nach dem Pivotieren die obigen
  129.            ; Gleichungen gelten.
  130.            (let ((r (/ (aref AA k l))))
  131.              (declare (type Zahl r))
  132.              ; Spalte l :
  133.              (progn
  134.                (dotimes (i m)
  135.                  (unless (eql i k)
  136.                    (setf (aref AA i l) (- (* r (aref AA i l))))
  137.              ) ) )
  138.              (setf (aref AC l) (- (* r (aref AC l))))
  139.              ; alles außer Zeile k, Spalte l :
  140.              (dotimes (j n)
  141.                (unless (eql j l)
  142.                  (let ((s (aref AA k j)))
  143.                    (dotimes (i m)
  144.                      (unless (eql i k)
  145.                        (setf (aref AA i j) (+ (aref AA i j) (* s (aref AA i l))))
  146.                    ) )
  147.                    (setf (aref AC j) (+ (aref AC j) (* s (aref AC l))))
  148.              ) ) )
  149.              (let ((s (aref AB k)))
  150.                (dotimes (i m)
  151.                  (unless (eql i k)
  152.                    (setf (aref AB i) (+ (aref AB i) (* s (aref AA i l))))
  153.                ) )
  154.                (setf AD (+ AD (* s (aref AC l))))
  155.              )
  156.              (when Extend
  157.                (locally (declare (type Zahl-Matrix AAE) (type Zahl-Vektor ACE))
  158.                  (dotimes (j m)
  159.                    (let ((s (aref AAE k j)))
  160.                      (dotimes (i m)
  161.                        (unless (eql i k)
  162.                          (setf (aref AAE i j) (+ (aref AAE i j) (* s (aref AA i l))))
  163.                      ) )
  164.                      (setf (aref ACE j) (+ (aref ACE j) (* s (aref AC l))))
  165.              ) ) ) )
  166.              ; Zeile k :
  167.              (progn
  168.                (dotimes (j n)
  169.                  (unless (eql j l)
  170.                    (setf (aref AA k j) (* (aref AA k j) r))
  171.                ) )
  172.                (setf (aref AB k) (* (aref AB k) r))
  173.              )
  174.              (when Extend
  175.                (locally (declare (type Zahl-Matrix AAE))
  176.                  (dotimes (j m)
  177.                    (setf (aref AAE k j) (* (aref AAE k j) r))
  178.              ) ) )
  179.              ; Element (k,l) :
  180.              (setf (aref AA k l) r)
  181.              ; Markierungen vertauschen:
  182.              (rotatef (aref Zeile k) (aref Spalte l))
  183.         )) )
  184.         ; Variablen Z nach unten schaffen, evtl. Matrix verkleinern:
  185.         (let ((elbar (make-array m :fill-pointer 0))
  186.               (not-elbar (make-array m :fill-pointer 0)))
  187.           ; elbar = Menge der eliminierbaren z,
  188.           ; not-elbar = Menge der nicht eliminierbaren z.
  189.           (dotimes (i m)
  190.             ; Suche Betragsmaximum in Zeile i:
  191.             (let ((hmax 0) (l nil))
  192.               (dotimes (j n)
  193.                 (when (eq (car (aref Spalte j)) 'XY)
  194.                   (let ((h (abs (aref AA i j))))
  195.                     (when (> h hmax) (setq hmax h l j))
  196.               ) ) )
  197.               (if l
  198.                 ; AA[i,l] war betragsmäßig maximal
  199.                 (progn
  200.                   (vector-push i not-elbar)
  201.                   (pivot i l)
  202.                 )
  203.                 ; triviale Zeile
  204.                 (if (zerop (aref AB i))
  205.                   ; Führe diese Zeile noch pro forma weiter; sie wird sich
  206.                   ; aber nicht weiter ändern, weil stets nur um XY-Spalten
  207.                   ; pivotiert werden wird und in diesen Spalten in Zeile i
  208.                   ; überall 0 steht.
  209.                   ; Nicht-Null-Elemente stehen nur in ZO-Spalten, und die
  210.                   ; verändern sich nicht mehr.
  211.                   (vector-push i elbar)
  212.                   ; Zeile fürs LP widerspruchsvoll -> unlösbar
  213.                   (return-from simplex NIL)
  214.           ) ) ) )
  215.           ; eliminierbare Zeilen in der ZZ-Matrix vormerken:
  216.           ; Die nicht eliminierbaren Zeilen haben mit den X ihre Plätze
  217.           ; getauscht, so daß also Spalte[(cdr Zeile[i])] = (ZO . i) gilt,
  218.           ; falls Zeile i nicht eliminierbar ( <==> (car Zeile[i]) = XY ).
  219.           (dotimes (i0h (fill-pointer elbar))
  220.             (let ((i0 (aref elbar i0h)))
  221.               (setf (aref Z i0) 0) ; nur wenn man Z[i0]=0 setzt, stimmen die anderen Z's!
  222.               (setf (aref ZZ i0 i0) T) ; z_i0 frei wählbar
  223.               (dotimes (ih (fill-pointer not-elbar))
  224.                 (let* ((i (aref not-elbar ih)) ; es muß (car Zeile[i])=XY sein
  225.                        (j (cdr (aref Zeile i)))) ; Spalte, mit der Zeile i pivotiert wurde
  226.                   (setf (aref ZZ i i0) (aref AA i0 j))
  227.           ) ) ) )
  228.           ; Zeilen streichen: (benutze dabei, daß alle not-elbar[i]>=i)
  229.           (dotimes (ih (fill-pointer not-elbar))
  230.             (let ((i (aref not-elbar ih)))
  231.               (unless (eql ih i)
  232.                 (dotimes (j n) (setf (aref AA ih j) (aref AA i j)))
  233.                 (setf (aref AB ih) (aref AB i))
  234.           ) ) )
  235.           (setq m (fill-pointer not-elbar)) ; neue Zeilenzahl = Anzahl der ZO-Spalten
  236.         )
  237.         ; Spalten umsortieren: XY links, ZO rechts bringen.
  238.         ; Damit werden am Schluß die nicht eliminierbaren z aus den X errechnet.
  239.         (let ((l 0) ; linke Spalte
  240.               (r (1- n))) ; rechte Spalte
  241.           (loop
  242.             (unless (< l r) (return))
  243.             (cond ((eq (car (aref Spalte l)) 'XY) (incf l))
  244.                   ((eq (car (aref Spalte r)) 'ZO) (decf r))
  245.                   (t ; Vertausche Spalten r und l
  246.                      (dotimes (i m) (rotatef (aref AA i l) (aref AA i r)))
  247.                      (rotatef (aref AC l) (aref AC r))
  248.                      (rotatef (aref Spalte l) (aref Spalte r))
  249.             )     )
  250.         ) )
  251.         ; Verstecke diese M Spalten vor dem Pivotieren:
  252.         (setq n (- n m))
  253.         ; Die Elemente AA[0..m-1,n..n+m-1], AC[n..n+m-1], Spalte[n,n+m-1]
  254.         ; werden erst am Schluß von Teil 6 wieder benötigt.
  255.         (let ((Zeile_save (copy-seq Zeile)))
  256.           (flet
  257.             ((SuchePivotZeile (AWZM l)
  258.                ; Bestimmt zur Auswahlmenge AWZM der zu untersuchenden Zeilen
  259.                ; und zur Spalte l die Zeile k so:
  260.                ; Vorausgesetzt ist, daß
  261.                ; für i∈AWZM gilt 0≤i<m und AB[i]≥0 und AA[i,l]>0.
  262.                ; k ist unter allen i∈AWZM dasjenige, für das der Quotient
  263.                ; AB[i]/AA[i,l] minimal ist. Sollte er =0 sein oder ist sowieso
  264.                ; Extend=true, so ist nachher Extend=true, und es wurde der Vek-
  265.                ; tor (AB[i]/AA[i,l], AAE[i,1]/AA[i,l], ..., AAE[i,m-1]/AA[i,l])
  266.                ; unter allen i∈AWZM lexikographisch minimiert. Falls AWZM leer
  267.                ; ist, ist NIL das Ergebnis.
  268.                (if (eql (fill-pointer AWZM) 0)
  269.                  NIL
  270.                  (let (k)
  271.                    (unless Extend ; versuche k passend zu wählen
  272.                      (let (hmax)
  273.                        (dotimes (ih (fill-pointer AWZM))
  274.                          (let* ((i (aref AWZM ih))
  275.                                 (h (/ (aref AB i) (aref AA i l))))
  276.                            (when (or (eql ih 0) (< h hmax)) (setq hmax h k i))
  277.                        ) )
  278.                        (when (zerop hmax)
  279.                          ; entarteter Fall
  280.                          (setq Extend T)
  281.                          (setq AAE (make-array `(,m ,m) :initial-element 0))
  282.                          (dotimes (i m) (setf (aref AAE i i) 1))
  283.                          (setq ACE (make-array m :initial-element 0))
  284.                    ) ) )
  285.                    (when Extend
  286.                      ; Der Entartungsfall liegt entweder seit längerem vor
  287.                      ; oder ist erst jetzt aufgetreten (dann war das alte k
  288.                      ; vielleicht schlecht).
  289.                      (let (hmax hmaxe)
  290.                        (dotimes (ih (fill-pointer AWZM))
  291.                          (let ((i (aref AWZM ih)))
  292.                            (let ((h (/ (aref AA i l)))
  293.                                  (he (make-array m)))
  294.                              (dotimes (j m) (setf (aref he j) (* (aref AAE i j) h)))
  295.                              (setq h (* (aref AB i) h))
  296.                              ; (h,he[1],...,he[M]) ist jetzt der Quotientenvektor
  297.                              (when (or (eql ih 0)
  298.                                        ; statt (< h hmax) jetzt lexikographisch
  299.                                        ; (h,he[1],...,he[M]) < (hmax,hmaxe[1],...,hmaxe[M])
  300.                                        ; vergleichen:
  301.                                        (or (< h hmax)
  302.                                            (and (= h hmax)
  303.                                                 (dotimes (ie m NIL)
  304.                                                   (let* ((he_ie (aref he ie))
  305.                                                          (hmaxe_ie (aref hmaxe ie)))
  306.                                                     (when (< he_ie hmaxe_ie) (return T))
  307.                                                     (when (> he_ie hmaxe_ie) (return NIL))
  308.                                    )   )   )    ) )
  309.                                (setq hmax h hmaxe he k i)
  310.                    ) ) ) ) ) )
  311.                    k
  312.             )) ) )
  313.             (let ((PZM (make-array m :fill-pointer 0)))
  314.               ; PZM = Menge der i mit AB[i]≥0
  315.               (loop ; Spalte AB mit positiven Zahlen füllen
  316.                 #|
  317.                 ; Rundungsfehler ausmerzen:
  318.                 (dotimes (ih (fill-pointer PZM))
  319.                   (let ((i (aref PZM ih)))
  320.                     (when (minusp (aref AB i)) (setq (aref AB i) 0))
  321.                 ) )
  322.                 |#
  323.                 (let ((NZM (make-array m :fill-pointer 0)))
  324.                   ; NZM = Menge der i mit AB[i]<0
  325.                   ; PZM und NZM neu berechnen:
  326.                   (let ((old-PZM-count (fill-pointer PZM))) ; alte Mächtigkeit von PZM
  327.                     (setf (fill-pointer PZM) 0)
  328.                     (dotimes (i m)
  329.                       (if (>= (aref AB i) 0)
  330.                         (vector-push i PZM)
  331.                         (vector-push i NZM)
  332.                     ) )
  333.                     ; Streiche die Zusatzspalten, wenn sich PZM echt vergrößert
  334.                     ; hat und die Entartung daher vielleicht verschwunden ist.
  335.                     (when (> (fill-pointer PZM) old-PZM-count) (setq Extend nil))
  336.                   )
  337.                   (when (eql (fill-pointer NZM) 0) (return)) ; alle AB[i]≥0 ?
  338.                   ; wähle p aus: p:= dasjenige i mit AB[i]<0,
  339.                   ; bei dem die Anzahl der AA[i,j]≥0 maximal ist.
  340.                   (let ((countmax -1) (p nil))
  341.                     (dotimes (ih (fill-pointer NZM))
  342.                       (let ((i (aref NZM ih))
  343.                             (count 0))
  344.                         (dotimes (j n) (when (>= (aref AA i j) 0) (incf count)))
  345.                         (when (> count countmax) (setq countmax count p i))
  346.                     ) )
  347.                     (when (eql countmax n)
  348.                       ; alle AA[p,j]≥0, aber AB[p]<0
  349.                       ; ==> Zeile fürs LP widerspruchsvoll ==> unlösbar
  350.                       (return-from simplex NIL)
  351.                     )
  352.                     ; Wähle l aus: maximales abs(AA[p,j]) unter allen j mit AA[p,j]<0
  353.                     ; (solche gibt es, da countmax<n).
  354.                     (let ((hmin 0) (l nil))
  355.                       (dotimes (j n)
  356.                         (let ((h (aref AA p j)))
  357.                           (when (< h hmin) (setq hmin h l j))
  358.                       ) )
  359.                       ; AWZM aufbauen:
  360.                       (let ((AWZM (make-array m :fill-pointer 0)))
  361.                         (dotimes (ih (fill-pointer PZM))
  362.                           (let ((i (aref PZM ih)))
  363.                             (when (> (aref AA i l) 0) (vector-push i AWZM))
  364.                         ) )
  365.                         (let ((k (SuchePivotZeile AWZM l)))
  366.                           (if (null k)
  367.                             ; Durch Pivotieren um AA[p,l] wächst PZM um
  368.                             ; mindestens das eine Element p, denn:
  369.                             ; für i∈PZM gilt AB[i]≥0 und wegen AZWM=ϕ
  370.                             ; auch AA[i,l]≤0, also nachher
  371.                             ; AB[i]=AB[i]-AB[p]*AA[i,l]/AA[p,l] ≥0,
  372.                             ;        ≥0    <0     ≤0     <0
  373.                             ; also wieder i∈PZM.
  374.                             ; Aber nachher auch AB[p]=AB[p]/AA[p,l] (<0/<0) >0,
  375.                             ; also p∈PZM.
  376.                             (progn
  377.                               (setq Extend nil) ; brauche die Zusatzspalten jetzt nicht mehr
  378.                               (pivot p l)
  379.                             )
  380.                             ; Pivotzeile k ausgewählt; fertig zum Pivotieren
  381.                             ; Durch Pivotieren um AA[k,l] wird PZM nicht
  382.                             ; verkleinert, denn:
  383.                             ; für i∈PZM gilt AB[i]≥0.
  384.                             ; Falls AA[i,l]≤0, gilt nachher
  385.                             ; AB[i]=AB[i]-AB[k]*AA[i,l]/AA[k,l] ≥0.
  386.                             ;        ≥0    ≥0     ≤0      >0
  387.                             ; Falls jedoch AA[i,l]>0, so gilt für i<>k nachher
  388.                             ; AB[i]=AA[i,l]*(AB[i]/AA[i,l]-AB[k]/AA[k,l])     ≥0
  389.                             ;          >0    ≥0 nach Wahl von k wegen i∈AWZM
  390.                             ; und für i=k nachher AB[k]=AB[k]/AA[k,l] (≥0/>0) ≥0.
  391.                             ; Also nachher stets wieder AB[i]≥0, d.h. i∈PZM.
  392.                             (pivot k l)
  393.                 ) ) ) ) ) )
  394.             ) )
  395.             ; Jetzt ist Extend=false, weil (fill-pointer PZM) sich gerade
  396.             ; auf m erhöht haben muß.
  397.             ; Ab hier sind alle AB[i]≥0, und diese Eigenschaft geht nicht verloren.
  398.             (loop
  399.               ; Suche ein l mit AC[l]<0 :
  400.               (let (l)
  401.                 (when
  402.                   (dotimes (j n T)
  403.                     (when (< (aref AC j) 0) (setq l j) (return NIL))
  404.                   )
  405.                   (return) ; alle AC[j]≥0 ==> lösbar
  406.                 )
  407.                 ; AWZM := Menge der i mit AA[i,l]>0 :
  408.                 (let ((AWZM (make-array m :fill-pointer 0)))
  409.                   (dotimes (i m)
  410.                     (when (> (aref AA i l) 0) (vector-push i AWZM))
  411.                   )
  412.                   (let ((k (SuchePivotZeile AWZM l)))
  413.                     (if (null k)
  414.                       ; alle AA[i,l]≤0 und AC[l]<0 ==> Spalte fürs DP widerspruchsvoll
  415.                       (return-from simplex NIL)
  416.                       (pivot k l)
  417.                       ; weiterhin alle AB[i]≥0.
  418.                       ; AD:=AD-AB[k]*AC[l]/AA[k,l] ≥ AD wird nicht verkleinert.
  419.                       ;         ≥0    <0     >0
  420.               ) ) ) )
  421.             )
  422.             ; Lösbar! Lösung bauen:
  423.             (let ((vollständig t))
  424.               (dotimes (i m)
  425.                 (let ((s (aref AB i))
  426.                       (index (cdr (aref Zeile i))))
  427.                   (setf (aref X index) 0 (aref Y index) s)
  428.                   (setq vollständig (and vollständig (> s 0)))
  429.               ) )
  430.               (dotimes (j n)
  431.                 (let ((s (aref AC j))
  432.                       (index (cdr (aref Spalte j))))
  433.                   (setf (aref X index) s (aref Y index) 0)
  434.                   (setq vollständig (and vollständig (> s 0)))
  435.               ) )
  436.               ; Die nicht eliminierbaren z-Werte berechnen sich aus dem
  437.               ; verstecken Teil von AA und AC und den Werten von X und Y :
  438.               (do ((j n (1+ j)))
  439.                   ((>= j (+ n m)))
  440.                 (let ((s (aref AC j)))
  441.                   (dotimes (i m)
  442.                     (setq s (+ s (* (aref AA i j) (aref X (cdr (aref Zeile_save i))))))
  443.                   )
  444.                   (setf (aref Z (cdr (aref Spalte j))) s)
  445.               ) )
  446.               (values T (list Y (- AD)) (list X Z ZZ AD) vollständig)
  447. ) ) ) ) ) ) )
  448.  
  449. (defun test-simplex (aufg)
  450.   (multiple-value-bind (m n) (values-list (array-dimensions (first aufg)))
  451.     (let ((x-list
  452.             (mapcar #'(lambda (i) (format nil "X[~D]" i))
  453.                     (make-sequence 'list n :initial-element 1 :update #'1+)
  454.           ) )
  455.           (y-list
  456.             (mapcar #'(lambda (i) (format nil "Y[~D]" i))
  457.                     (make-sequence 'list n :initial-element 1 :update #'1+)
  458.           ) )
  459.           (z-list
  460.             (mapcar #'(lambda (i) (format nil "Z[~D]" i))
  461.                     (make-sequence 'list m :initial-element 1 :update #'1+)
  462.          )) )
  463.       ; Aufgabe ausgeben:
  464.       (multiple-value-bind (A b c) (values-list aufg)
  465.         (format t "~%~%Lineares Problem:")
  466.         (dotimes (i m)
  467.           (let ((zeile (make-array n)))
  468.             (dotimes (j n) (setf (aref zeile j) (aref A i j)))
  469.             (format t "~%~{~1{~S * ~A~:}~^ + ~} = ~S"
  470.                       (remove 0 (mapcar #'list (coerce zeile 'list) y-list)
  471.                                 :key #'first :test #'=
  472.                       )
  473.                       (aref b i)
  474.         ) ) )
  475.         (format t "~%mit ~{~A~^, ~} ≥ 0" y-list)
  476.         (format t "~%Minimiere ~{~1{~S * ~A~:}~^ + ~}."
  477.                   (remove 0 (mapcar #'list (coerce c 'list) y-list)
  478.                             :key #'first :test #'=
  479.         )         )
  480.         (format t "~%Dualproblem:")
  481.         (dotimes (j n)
  482.           (let ((spalte (make-array m)))
  483.             (dotimes (i m) (setf (aref spalte i) (aref A i j)))
  484.             (format t "~%~A = ~:{~S * ~A + ~} ~S ≥ 0"
  485.                       (elt x-list j)
  486.                       (remove 0 (mapcar #'list (coerce spalte 'list) z-list)
  487.                                 :key #'first :test #'=
  488.                       )
  489.                       (aref c j)
  490.         ) ) )
  491.         (format t "~%Minimiere ~{~1{~S * ~A~:}~^ + ~}."
  492.                   (remove 0 (mapcar #'list (coerce b 'list) z-list)
  493.                             :key #'first :test #'=
  494.         )         )
  495.       )
  496.       (let ((lösung (multiple-value-list (apply #'simplex aufg))))
  497.         ; Lösung ausgeben:
  498.         (if (first lösung)
  499.           (progn
  500.             (format t "~%~%~:[Eine~;Die einzige~] Lösung ist:" (fourth lösung))
  501.             (let ((LP-lösung (second lösung)))
  502.               (format t "~%~{~1{~A = ~S~:}~^, ~}~%Minimum = ~S"
  503.                         (mapcar #'list y-list (coerce (first LP-lösung) 'list))
  504.                         (second LP-lösung)
  505.             ) )
  506.             (let ((DP-lösung (third lösung)))
  507.               (let ((Z (second DP-lösung))
  508.                     (ZZ (third DP-lösung)))
  509.                 (dotimes (i m)
  510.                   (format t "~%~A = " (elt z-list i))
  511.                   (if (aref ZZ i i)
  512.                     (format t "beliebig")
  513.                     (format t "~S~:{ + ~S * ~A~}"
  514.                               (aref Z i)
  515.                               (let ((L nil))
  516.                                 (dotimes (j m)
  517.                                   (when (aref ZZ j j)
  518.                                     (unless (zerop (aref ZZ i j))
  519.                                       (push (list (aref ZZ i j) (elt z-list j))
  520.                                             L
  521.                                 ) ) ) )
  522.                                 (nreverse L)
  523.                     )         )
  524.               ) ) )
  525.               (format t "~%~{~1{~A = ~S~:}~^, ~}~%Minimum = ~S"
  526.                         (mapcar #'list x-list (coerce (first DP-lösung) 'list))
  527.                         (fourth DP-lösung)
  528.           ) ) )
  529.           (format t "~%unlösbar")
  530.         )
  531.         lösung
  532. ) ) ) )
  533.  
  534. (defun test ()
  535.   (mapcar #'test-simplex
  536.     '((#2A((1 0 1 -1  -1  3)
  537.            (0 1 2 -1 -1/2 1))
  538.        #(0 0)
  539.        #(0 0 -1 -1 -3 8)
  540.       )
  541.       (#2A((1 -1  -1  3 1 0)
  542.            (2 -1 -1/2 1 0 1))
  543.        #(0 0)
  544.        #(-1 -1 -3 8 0 0)
  545.       )
  546.       (#2A((1 1 0 1)
  547.            (2 0 1 0))
  548.        #(1 1)
  549.        #(-1 1 -2 0)
  550.       )
  551.       (#2A((1  1 2 1 3 -1  0)
  552.            (2 -2 3 1 1  0 -1))
  553.        #(4 3)
  554.        #(2 3 5 2 3 0 0)
  555.       )
  556.       (#2A(( 1  8  0  -1  0  0  0)
  557.            (-3 -7 -20  0 -1  0  0)
  558.            ( 1  0  0   0  0 -1  0)
  559.            ( 0 -1  0   0  0  0 -1)
  560.            ( 1  1  1   0  0  0  0))
  561.        #(3 -6 2/5 -1/2 1)
  562.        #(5 2 1/4 0 0 0 0)
  563.       )
  564.       (#2A((1 2  1  1 -2)
  565.            (1 1 -4 -2 -3)
  566.            (1 2  5 -4  6))
  567.        #(4 2 3)
  568.        #(3 5 4 5 6)
  569.       )
  570.       (#2A((1  2 -1  1)
  571.            (2 -2  3  3)
  572.            (1 -1  2 -1))
  573.        #(0 9 6)
  574.        #(-3 1 3 -1)
  575.      ))
  576. ) )
  577.  
  578.