Rozwiązywanie zadań programowania liniowego z wbudowaną strukturą sieciowo-transportową

ROCZNIKI POLSKIEGO TOWARZYSTWA MATEMATYCZNEGO Seria III: MATEMATYKA STOSOWANA XXXIII (1991)

K

Z

Warszawa

Rozwiązywanie zadań programowania liniowego z wbudowaną strukturą sieciowo-transportową

( Praca wpłynęła do Redakcji 7.2.1989)

S treszczen ie. W pracy przedstawiony jest uniwersalny algorytm rozwiązywania za- gadnień optymalnej dystrybucji w sieci transportowej, poddanej dodatkowym ogranicze- niom liniowym, czyli tzw. zadań programowania liniowego z wbudowaną strukturą sie- ciową. Algorytm ten funkcjonuje na zasadzie pierwotnej metody sympleksowej i opiera się na dekompozycji bazy na cztery bloki. Czynnikiem decydującym o efektywności algorytmu jest sposób realizacji operacji z udziałem bloku sieciowego. Dlatego szczególny nacisk po- łożony jest w pracy na zaprojektowaniu struktur danych uwzględniających specyfikę tego bloku i pozwalających wykorzystać ją w pełni na poziomie implementacji.

1. Wstęp. Wśród modeli programowania liniowego najliczniejsze zasto- sowania praktyczne znalazły liniowe modele transportowe, lub szerzej rozu- miane sieciowe modele dustrybucji. Każdy z tych modeli charakteryzuje się tym, że opisuje sieć transportową, której obrazem jest pewien graf zoriento- wany. W lukach sieci płynie, zgodnie z ich orientacją, jednorodny produkt w ilościach nie przekraczających (w jednostce czasu) danych przepustowości łuków. Produkt ten jest emitowany i pochłaniany przez węzły. Jeśli dany węzeł emituje w jednostce czasu e jednostek produktu a pochłania p, to z punktu widzenia jego funkcjonowania istotny jest tylko bilans w postaci różnicy e — p. Tak więc każdy węzeł sieci scharakteryzowany jest przez dany bilans.

Przez dopuszczalny schemat dystrybucji rozumie się dowolny układ prze-

pływów łukowych spełniających warunki bilansowe i ograniczenia przepu-

stowości łuków. Jeśli każdemu łukowi przyporządkowany jest współczynnik

Praca była wykonana w ramach programu resortowego RP.I.02 - “Teoria Sterowania

i Optymalizacji Ciągłych Układów Dynamicznych i Procesów Dyskretnych ”

liczbowy zwany umownie kosztem jednostkowym przepływu, to każdy sche- mat dystrybucji można ocenić za pomocą globalnego kosztu dystrybucji, będącego sumą łukowych liniowych kosztów przepływu. Tak więc najlepszy schemat dystrybucji może być utożsamiany ze schematem odpowiadającym minimalnemu globalnemu kosztowi dystrybucji. Aby go znaleźć, wystarczy rozwiązać odpowiednie zadanie programowania liniowego polegające na mi- nimalizacji kosztu globalnego przy warunkach bilansowych i ograniczeniach przepustowościowy ch.

Z uwagi na specyficzną postać macierzy (macierz całkowicie unimodu- larna) zadania optymalnej dystrybucji w sieci należą do podklasy zadań pro- gramowania liniowego dysponującej własną rozbudowaną teorią i własnymi metodami wykorzystującymi charakterystyczne cechy zadań tej klasy. W y- korzystanie tych cech znacznie ułatwia (zwłaszcza pod względem numerycz- nym) rozwiązywanie i zwiększa efektywność algorytmów (por. Grigoriadis [5], Jabłonowski, Zorychta [7]).

Zagadnienia praktyczne rzadko jednak przystają do modeli transporto- wych i dystybucyjnych w postaci czystej. W modelach praktycznych oprócz warunków bilansowych i przepustowościowych pojawiają się zwykle dodat- kowo ogólniejsze ograniczenia. Dla przykładu, w zagadnieniach transporto- wo-lokalizacyjnych będą to warunki całowitoliczbowości zmiennych wskazu- jących lokalizacje, lub zmienne (zależne od lokalizacji) górne ograniczenia na przepływy łukowe. W zagadnieniach wielokryterialnych do równań bi- lanswych dołączone będą przeważnie równania, w które uwikłane są kolejne funkcje celu rozpatrywane w modelu. Mogą to być również inne zależności liniowe, np.-obejmujące przepływy łukowe z pewnych wyróżnionych rejonów sieci w celu ograniczenia obciążenia transportowego itp.

Takie rozszerzenie modelu dystrybucji w sieci prowadzi do zadań pro- gramowania liniowego z tzw. wbudowaną strukturą sieciową. Jeśli struktura sieciowa stanowi tyko niewielką część całego modelu, to nie ma powodu trak- tować takie zadanie inaczej niż ogólne zadanie programowania liniowego.

Jeżeli jednak wbudowana część sieciowa stanowi zasadniczą część modelu i pozostałe warunki są nieliczne w stosunku do liczby równań bilansowych, to warto zastanowić się w jaki sposób wykorzystać tę specyficzną struk- turę zadania w celu zwiększenia efektywności algorytmu i poprawienia jego własności numerycznych.

Celem tej pracy jest przedstawienie uniwersalnego algorytmu rozwiązy-

wania tego typu zadań (algorytm SON, por. Glover, Klingman [3,4]). Algo-

rytm ten oparty jest na dekompozycji bieżącej bazy na cztery bloki, z któ-

rych blok odpowiadający równaniom i zmiennym sieciowym zdecydowannie

dominuje rozmiarami nad pozostałymi blokami. Cała numeryka tego algo-

rytmu związana jest z tzw. jądrem, czyli blokiem bazy odwrotnej odpowia-

dającym (w pewnym uproszczeniu) zmiennym i równaniom niesieciowym.

Rozwiązywanie zadań programowania liniowego 13 Różne sposoby prowadzenia jądra lub jego reprezentacji, odnawiania go przy zmianie bazy oraz realizacji operacji z jego udziałem mogą prowadzić do istotnie różnych wersji algorytmu. Wersje te różnić się będą między sobą wła- snościami numerycznymi oraz kosztem realizacji poszczególnych operacji.

W niniejszej pracy zaproponowany jest sposób przekształcania i odnawia- nia jądra danego explicite, oparty na metodzie eliminacji Gaussa-Jordana.

Wprawdzie sposób ten nie gwarantuje stabilności prowadzonych obliczeń, jednak niewielkie wymiary jąder dla jakich jest przewidziany oraz możli- wość stosowania podwójnej precyzji w obliczeniach numerycznych stwarzają nadzieję, że prezentowany algorytm będzie użyteczny w rozwiązywaniu dość obszernej klasy zadań praktycznych. Niewątpliwą jego zaletą jest prostota realizacji, co szczególnie preferuje go przy implementowaniu algorytmu na mikrokomputerach typu IBM P C /X T lub AT, kiedy to zarówno czas obli- czeń jak i dostępna dla obliczeń pamięć muszą być brane pod uwagę (por.

Ogryczak, Studziński, Zorychta [8]).

2. Zadanie program ow ania liniowego z w budow aną strukturą sieciow ą. Rozważmy zadanie programowania liniowego postaci

(2.1) min {err : A x = b , 0 < x < g}

w którym wszystkie dane: A = («ij)> b = (6;), c = (cj), g = (gj) (i E 1, A:, j E .1,/) oraz zmienne decyzyjne x = (xj) (j E 1,/) są rzeczywiste, przy czym gj < +oo.

Załóżmy, że macierz A dzieli się na następujące bloki

( 2 ^. 2 ^{) A} _A

^A

l n

A

Odpowiednio do tego podziału dzielą się również na części wektory 6, c, g oraz x. Podział ten wraz z wymiarami bloków przedstawiony jest schema- tycznie na rys. 1.

<--- 1 ---►

n--- k --- q

Xn X L

Cn CL

An n An l bN

Al n Al l h

£ N CL

Rys. 1

Zadanie (2.1) jest zadaniem programowania liniowego z wbudowaną

strukturą sieciową jeśli jego macierz można podzielić na bloki w taki sposób,

by blok A

miał strukturę sieciową.

Będziemy mówili, że dana macierz H ma strukturę sieciową jeśli spełnia następujące warunki:

(i) każdy jej element jest równy 0, 1 lub -1;

(ii) w każdej kolumnie ma nie mniej niż jeden i nie więcej niż dwa ele- menty różne od zera;

(iii) jeśli w kolumnie są dwa elementy różne od zera, to są one różnych znaków.

Danej macierzy H o strukturze sieciowej można przyporządkować sieć transportową G w następujący sposób. Każdemu wierszowi macierzy od- powiada wzajemnie jednoznacznie jeden węzeł sieci. Dlatego wiersze takiej macierzy będziemy utożsamiali z węzłami odpowiadającej jej sieci. Kolumny macierzy H będą natomiast identyfikowane z lukami. Jeśli kolumna zawiera dwa elementy niezerowe ( — 1 i +1), to odpowiadający jej luk zaczyna się w węźle wskazanym przez —l a kończy w węźle wskazanym przez +1. Takie luki nazywać będziemy lukami zwyczajnymi. Jeśli kolumna zawiera dokład- nie jeden element różny od zera, wskazujący jeden z końców luku ( —1 - początek, a +1 - koniec luku), to jako drugi koniec luku traktowany jest pewien węzeł spoza sieci. Takie luki nazywać będziemy lukami uzupełniają- cymi.

Macierz zadania z wbudowaną strukturą sieciową zawiera zatem blok A

odpowiadający czystemu problemowi dystrybucji w sieci transporto- wej, gdyż z dokładnością do kolumn zawierających po jednym elemencie niezerowym macierz A

jest macierzą incydencji węzłów i łuków odpo- wiadającej jej sieci. Następujący lemat wyraża znany fakt dający pewną charakterystykę macierzy o strukturze sieciowej (por. Ilu [6]):

L

1. Każda macierz o strukturze sieciowej jest całkowicie unimo- dularna , tzn. każda jej podmacierz kwadratowa ma wyznacznik równy 0, 1 lub —1.

Fakt ten ma duże znaczenie dla rozwiązywania zadań programowania liniowego z wbudowaną strukturą sieciową. Dotyczy to zarówno budowy al- gorytmów jak i ich własności numerycznych. Pełne wykorzystanie właściwo- ści struktury sieciowej w istniejących algorytmach programowania liniowego pozwala zmniejszyć ich koszt, a także w znacznej mierze eliminuje wpływ części sieciowej na produkcję błędów numerycznych.

3. P ierw otn a m etod a sym pleksow a. Dla rozwiązania zadania (2.1)

wykorzystany zostanie wariant pierwotnej metody sympleksowej znany pod

nazwą SUB (por.Zorychta, Ogryczak [10]) i dostosowany do zadań z górnymi

ograniczeniami. W tym celu założymy, że rząd macierzy A równy jest licz-

bie jej wierszy (rankA = k ). Podstawę tej metody stanowi zmieniająca się z

Rozwiązywanie zadań programowania liniowego 15 iteracji na iterację baza B , reprezentująca bieżące bazowe rozwiązanie dopu- szczalne. Baza ta nie definiuje jednak rozwiązania bazowego jednoznacznie.

Dla jego określenia niezbędne jest jeszcze zdefiniowanie wartości zmiennych niebazowych, z których każda jest równa bądź swemu dolnemu (0), bądź gór- nemu ( gj ) ograniczeniu. Tak więc ta sama baza może reprezentować pewną liczbę rozwiązań bazowych, zróżnicowanych w wyniku przypisywania zmien- nym niebazowym różnych układów wartości ze zbioru ich dolnych i górnych ograniczeń.

Kolejne kroki metody SUB można krótko opisać w sposób następujący.

0. Inicjalizacja. Zakładamy, że dane jest początkowe bazowe rozwiązanie dopuszczalne x = (I

^

)? gdzie

(3.1) afo — B _1(6 — N

T

)? 0 < afu < gu

a I

* (* G 1,1 —k) mają przypisane w pewien sposób wartości swych ograniczeń 0 lub </

.

Podziały wektorów x = (^

,^

), fj = (<7

,#

)>

= (

,

) odpo- wiadają podziałowi macierzy A zadania (2.1) na część bazową i niebazową A = (B ,N ).

1. Wycena. Znajdujemy wektor dualny

(3.2) w =

B -1

oraz wartości kosztów zredukowanych zmiennych niebazowych

(3.3) ć =

— i N

2. Test optymalności i wybór zmiennej wchodzącej do bazy.

(a) Jeżeli

(3.4)

> 0 dla = 0

< 0 dla = gNj j G 1, / — k to bieżące rozwiązanie bazowe jest optymalne; STOP.

(b) W przeciwnym wypadku wybieramy zmienną niebazową ,

speł- niającą warunek

(3.5) (

< 0 i JF

= 0) lub (

> 0 i £

= </Ns) jako zmienną wchodzącą do bazy.

3. Reprezentacja bazowa kolumny wchodzącej do bazy. Wyznaczamy re- prezentację a kolumny o wchodzącej do bazy (odpowiadającej wybranej zmiennej

«)

(3.6) o = B 1 o

(3.9) a = f 1 jeśli xNs =

\ - 1 jeśli xNs =

Ą. Test nieograniczoności i zmiana rozwiązania bazowego. Wyznaczamy maksymalne A spełniające warunki

(3.7) 0 < A < gNs

0 < x 'D i< g Bi , i e l , k gdzie

(3.8) — *EB* CrActBi i *^Ns — T o'A są nowymi wartościami zmiennych decyzyjnych oraz

= 0 9 ns

(a) Jeśli A = -foo, to funkcja celu jest nieograniczona z dołu i zadanie nie ma rozwiązania optymalnego; STOP.

(b) Jeśli A = gNs < +oo, to zmiennym bazowym oraz zmiennej nieba- zowej xns przypisujemy nowe wartości zgodnie z wzorami (3.8) i bez zmiany bazy przechodzimy do punktu 2.

(c) W pozostałych przypadkach zmieniamy wartości zmiennych decy- zyjnych zgodnie z wzorami (3.8) i przechodzimy do punktu 5.

5. Zmiana bazy i odnowienie reprezentacji jej odwrotności. Do bazy wcho- dzi kolumna odpowiadająca zmiennej niebazowej x n s

miejsce kolumny bazowej odpowiadającej zmiennej x s r , która w wyniku zmiany rozwiąza- nia w punkcie 4 osiągnęła swoje dolne (0) lub górne (</Br) ograniczenie, i w wyniku tego przesunięta została do odpowiedniej grupy zmiennych niebazo- wych. Odwrotność nowej bazy B ' wyrazi się następująco

(3.10) (B ')_1 = P B -1

gdzie P oznacza transformację eliminacyjną Gaussa-Jordana, tzn.

P * (^4, . • • y ^r — 1 j ^r+1 ? • • • i &k) — I

(ej oznacza i-ty wektor jednostkowy, I - macierz jednostkową.)

Powyższy opis prezentuje w skrócie jedynie zasadę działania metody SUB. Pewne narzędzia służące w miarę efektywnej implementacji tej metody dla zadań z wbudowaną strukturą sieciową przedstawione będą w kolejnych punktach pracy.

4. P rzestaw ien ia w ierszy i kolum n przy operacjach m acierzo- w ych . Analiza opisanego w poprzednim punkcie algorytmu wskazuje, że operacjami decydującymi o koszcie jednej jego iteracji są operacje wyszcze- gólnione w punktach 1 i 3 opisu, z udziałem macierzy B -1 . Są to operacje postaci

(4.1) y = B - ' f

Rozwiązywanie zadań programowania liniowego 17

w = e B 1

Operacje te równoważne są rozwiązaniu odpowiadających im układów rów- nań

(4.2) B y = f

tu — e

Dlatego, w zależności od potrzeb będziemy posługiwali się jedną lub drugą formą określania uwikłanych w te związki niewiadomych y i w.

Mechanizmy, jakie będą w tym celu wykorzystywane niekoniecznie będą odwoływać się do macierzy B -1 lub B, lecz będą posługiwać się macierzami, które różnią się od nich jedynie pewnymi permutacjami wierszy i kolumn. W związku z tym ustalimy najpierw pewne standardowe, charakterystyczne dla danej iteracji uporządkowania wierszy i kolumn B, a następnie prześledzimy co się dzieje i jak to wpływa na wyznaczanie y i w ze związków (4.1) lub (4.2) gdy dopuścimy pewne zmiany tych uporządkowań.

Podstawową charakterystykę bieżącej iteracji sympleksowej stanowić bę- dzie uporządkowane pole numerów kolumn bazowych

(4.3) NrB = [B l,B '2,...,B k]

oraz funkcja charakterystyczna zmiennych nieba.zowych ( . .v i

o f 0 jeśli j 4 NrB oraz Xj = 0 (4.4) Chi[j] = | G jeśli j \ NrB oraz aJj = Qj

Indeks Bi wskazuje, że kolumna Asi macierzy A znajduje się w bazie na i-tej pozycji, a to oznacza, że Ti-1 A Bi = £i (i G l,fc), gdzie ej oznacza i-ty wektor jednostkowy.

Jako standardowe uporządkowanie równań, a więc wierszy macierzy A , a w konsekwencji i wierszy jej podmacierzy B i N przyjmujemy dla każdej iteracji uporządkowanie naturalne [1 ,2 ,..., Ar]. Standardowe uporządkowa- nie kolumn macierzy B będzie w każdej iteracji zdefiniowane przez aktualne pole numerów kolumn bazowych NrB.

Zanim scharakteryzujemy wpływ jaki ma na wykonanie operacji (4.1) - (4.2) zmiana kolejności wierszy lub kolumn macierzy B rozważymy pewne uogólnienie tych związków w postaci równania macierzowego D E = F. Rów- naniu temu równoważne jest równanie P D Q Q _1E R = P F R , gdzie P, Q, R są macierzami permutacyjnymi odpowiednich wymiarów.

W naszym algorytmie najczęściej spotykane przypadki tego związku są następujące:

(1°) P B Q Q ^ B - 1? - 1 = I

(2°) P B Q Q -1 a,* = P /

(3°) w Q Q -1 B R = eR

Warto odnotować następujące, wynikające z rozważenia tych przypadków obserwacje:

ad 1° Kolumnom macierzy B odpowiadają wzajemnie jednoznacznie wiersze macierzy B -1 i zmiana kolejności jednych z nich pociąga za sobą automatyczną zmianę kolejności drugich. W szczególności, jeśli B ma stan- dardowo uporządkowane wiersze w kolejności [1,2,•••,&] oraz kolumny w kolejności [J91, B2 , . . . , Bk], to wiersze macierzy B -1 odpowiadają kolum- nom B w kolejności [BI, B 2 ,..., Bk] a kolumny równaniom układu w ko- lejności [1, 2, . . . , k].

ad 2° Jeśli rozwiązujemy układ Ba* = / , to zmianie kolejności wierszy B musi towarzyszyć taka sama zmiana kolejności prawych stron układu, a zmiana kolejności kolumn B jest ściśle związana z analogiczną zmianą kolejności odpowiadających im niewiadomych.

ad 3° Jeśli rozwiązujemy układ wB = e, to zmianie kolejności kolumn B musi towarzyszyć odpowiednia zmiana kolejności współrzędnych wektora e, a zmiana kolejności wierszy B musi być związana z analogiczną zmianą kolejności niewiadomych w i.

Spostrzeżenia te stanowić będą podstawę dla budowy słowników i wskaź- ników pozwalających realizować poprawnie omawiane operacje mimo niekrę- pującego podejścia do kolejności wierszy i kolumn uczestniczących w tych operacjach macierzy.

5. D e k o m p o zy cja op era cji bazow ych Jak to już zostało stwierdzone, szczególna struktura zadania (2.1) znajduje swoje odbicie w postaci macierzy A , która w naturalny sposób dzieli się na cztery bloki (patrz. (2.2)). Podział ten można również odnieść do podmacierzy macierzy A , w szczególności do dowolnej bazy

W związku z tym dla zrealizowania operacji (4.1) (lub (4.2)) z udziałem macierzy B -1 (lub B ) posłużymy się wzorami na odrotność macierzy blo- kowej. Z uwagi na to, że macierz B ^ n nie musi być nieosobliwa (ani nawet kwadratowa), założymy, że z dokładnością do kolejności wierszy i kolumn macierz B jest podzielona na bloki

w taki sposób, że B u jest kwadratową, nieosobliwą podmacierzą macierzy B n h , maksymalnego stopnia (stopień B u jest równy rzędowi B^yv) • (5.1)

(5.2)

Rozwiązywanie zadań programowania liniowego 19 Odwrotność tak podzielonej na bloki macierzy B określona jest nastę- pująco

(5.3) B " 1 = BH1 + B ^ B u V ^ B ^ B l i 1 - B n 'B u V " 1 - V - 'B u B f , 1 V " 1

gdzie V = B 22 — Ta postać B 1 nie będzie wykorzystywana explicite, a iedynie posłuży nam do wyprowadzenia wzorów na rozwiązania układów (4.1) lub (4.2).

Wyznaczanie wektora cen dualnych. Przyjmijmy, że w = (wi,W2), przy czym podział na części w\ i w 2 odpowiada podziałowi wierszy B na część górną (odpowiadającą B u ) i dolną (odpowiadającą B 21), z zachowaniem odpowiedniej do kolejności wierszy macierzy kolejności współrzędnych w .

Część bazowa c b wektora c podzieli się również na dwie części c b = (

,

2)> jednak ten podział i kolejność współrzędnych odpowiadają po- działowi kolumn macierzy B na część sieciową B u i niesieciową B 12.

Korzystając z drugiego z wzorów (4.1) uzyskamy następujące wzory na obie części wektora w

(5.4) w2 = ( cb 2 ~

B ^ B ^ V " 1 wi = icm — W2B 21 jB^j1

Gdybyśmy dysponowali prostym i szybkim algorytmem obliczania ilo- czynu eB -1 oraz gdyby była dana explicite macierz V -1 , to wektor w mó- głby być obliczony według następującego schematu

(5.5) 1. 5 = CBl&ii

2. ^t = ^{CB2 ~} £'Bi2 3. u>2 = /V -1 4. s = cb i — W2B 21 5. Wi = -sBj"/

gdzie s i t są polami roboczymi o odpowiedniej długości. Narzędzia dla reali- zacji tych formuł przedstawione będą w kolejnych punktach pracy. Zwróćmy uwagę, że jeśli rozszerzymy macierz układu o wiersz współczynników funkcji celu, tzn. jeśli przyjmiemy

B 22 - ^{B 22 o} _{CB2 1} ^{B 21 -}

oraz V = B 22 — B 2i-B111B i2, to

(5.6) D = V -1 =

B21

c b i

V " 1 0

— Wo 1

B12 — (B i 2 , 0 )

Tym samym, mając w każdej iteracji macierz D , mamy jednocześnie V " 1 i część W2 wektora dualnego. Aby wyznaczyć część pierwszą W\ wystarczy zrealizować punkty 4 i 5 schematu (5.5).

Wyznaczanie reprezentacji bazowej kolumny wchodzącej do bazy. Zał- óżmy, że a = (0 1 , 0:2) jest kolumną macierzy A odpowiadającą zmiennej xs wchodzącej do bazy. Podział na części 01 i 02 odpowiada podziałowi zbioru wierszy B na część górną, odpowiadającą B u , i część dolną, odpo- wiadającą B 2 1, z zachowaniem odpowiedniej do kolejności wierszy kolejności współrzędnych o. Wynikiem mnożenia B -1 o będzie wektor o = (0 5 1 , ^ 2) podzielony na części odpowiadające podziałowi i kolejności kolumn w ma- cierzy B. Korzystając z piewszego z wzorów (4.1) uzyskamy następujące wyrażenia na oć b \ i 0 5 2 ^

(5.7) a B2 = V _1(o2 - B2iB{11ai)

«

= ^ ^ ( a i — B12ZFB2)

Bezpośredni dostęp do macierzy V -1 oraz prosty i szybki algorytm realizu- jący mnożenia typu B ^ 1/ ułatwiłby znacznie postępowanie prowadzące do wyznaczenia a. Postępowanie to mogłoby przebiegać według następującego schematu:

1. s — _{B n V} 2. t = a 2 - B 2

3. &B2 = ^{> 7}

4.

=

- B 12

5.

B r /*

Struktury i algorytmy umożliwiające realizację tych formuł przedstawione będą w kolejnych punktach opracowania.

Algorytmy (5.5) i (5.8) mogą być stosowane nie tylko do wyznaczania cen dualnych i rozkładów kolumn wprowadzanych do bazy, lecz również do odtwarzania całych (lub w części) wierszy i kolumn B-1 . Jest to niezbędne przy odnawianiu odwrotności bazy kiedy to dla odnowienia jądra trzeba znać wiersz główny transformacji eliminacyjnej (por. Zorychta, Ogryczak [10]) oraz wiersz i kolumnę, o które ewntualnie ma być uzupełnione jądro.

Formuły umożliwiające wyznaczanie wierszy i kolumn B -1 są następujące:

(5.9) Pi = (B -1 ); = e;B -1

Pj = (B -1 )J = B “ 'ej

gdzie indeks u dołu oznacza wiersz, u góry - kolumnę macierzy B _1, a e;, ej są wektorami jednostkowymi.

6. Bazowe drzewo rozpinające. Niech G° będzie grafem nieskiero-

Rozwiązywanie zadań programowania liniowego 21 wanym, który powstał z grafu odpowiadającego sieci transportowej przez zaniedbanie orientacji jego łuków. Kolejny lemat przedstawia pewien fakt, który stanowić będzie podstawę operacji z macierzą B _1.

L

2. Podzbiór kolumn macierzy A jv jv jest liniowo niezależny wtedy i tylko wtedy, jeśli zbiór odpowiadająych mu w grafie G° krawędzi nie zawiera cyklu.

Pozwoli to na zastąpienie operacji z udziałem B ^1 odpowiednimi ope- racjami na acyklicznym podgrafie grafu G'°, reprezentującym bazę. Będzie to tzw. bazowe drzewo rozpinające (lub po prostu drzewo). Oznaczmy go znakiem T. Zbiór węzłów tego drzewa składać się będzie z wszystkich wę- złów grafu G. Do węzłów tych zaliczymy ponadto jeden sztuczny węzeł, tzw.

korzeń główny. Ponieważ założyliśmy, że A

składa się z m wierszy, więc bazowe drzewo rozpinające składać się będzie z m + 1 węzłów { 0 , 1, . . . , m }, gdzie 0 reprezentuje korzeń główny, oraz z m krawędzi.

Węzły drzewa odpowiadające wierszom przynależnym do dolnej części bazy (B 21, B 22) nazywać będziemy korzeniami pomocniczymi. Każdy ko- rzeń pomocniczy jest podłączony do korzenia głównego za pomocą krawędzi pomocniczej, którą oznaczać będziemy symbolem EA.

Do drzewa włączamy wszystkie krawędzie odpowiadające łukom zwy- czajnym (czyli kolumnom sieciowym o dwóch elementach różnych od zera), należącym do części B u bazy. Część B u na ogół zawiera również łuki uzu- pełniające (kolumny o jednym elemencie różnym od zera). Łuki uzupełnia- jące mogą być dwóch rodzajów: (1) łuki, które są lukami uzupełniającymi w A

N , oraz (2 ) łuki zwyczajne w Ajvjv> dla których w wyniku podziału ma- cierzy na bloki jeden element niezerowy znalazł się w górnej a drugi w dolnej części macierzy B. Łuk uzupełniający pierwszego rodzaju reprezentowany jest w drzewie przez krawędź łączącą korzeń główny z węzłem wskazanym przez element niezerowy w odpowiadającej kolumnie. Łuk uzupełniający drugiego rodzaju reprezentuje krawędź łączącą oba węzły wskazane przez elementy niezerowe w kolumnie i osobliwość tej sytuacji polega jedynie na tym, że jeden z końców tej krawędzi jest korzeniem pomocniczym.

Pojawiają się teraz trzy pytania, na które postaramy się odpowiedzieć:

1. Jakiej struktury użyć dla reprezentacji bazowego drzewa rozpinają- cego, by była ona wygodna dla rozwiązywania układów

(6.1) B n x = f

w B n = e

2. Jakie są algorytmy rozwiązywania tych układów i jaki jest ich koszt ? 3. Jak odnawiać strukturę drzewa gdy zmienia się ono o jeden łuk w wyniku zmiany bazy ?

Opis struktury drzewa. Dla opisu drzewa wykorzystamy układ tablic

(m + l)-elementowych postaci Tabl [0 : m], z których jedne opisują jed- noznacznie strukturę drzewa, inne umożliwiają szybki przegląd jego wierz- chołków, a jeszcze inne ułatwiają rozwiązywanie układów (6.1) i odnawianie drzewa przy zmianie bazy.

Drzewo T jest drzewem rozpinającym w grafie G°, a więc orientacja łuków obowiązująca w grafie G — w drzewie T nie obowiązuje. Będzie ono natomiast miało swoją własną, naturalną orientację drzewa, tzn. jeśli do drzewa należy krawędź ( k , /) oraz węzeł / jest bardziej odległy od korzenia niż węzeł k , to węzeł k jest poprzednikiem węzła / (a węzeł / następnikiem węzła k) w drzewie T. Przez odległość danego węzła od korzenia rozumiemy liczbę krawędzi tworzących drogę łączącą ten węzeł z korzeniem.

Pierwszą tablicą opisującą budowę samego drzewa jest tablica poprzed- ników:

Pred[ 0]d=f —1

Pred[i\ = poprzednik węzła i w drzewie T, i = 1,. . . , m

Każdy węzeł k drzewa (z wyjątkiem korzenia) odpowiada jednoznacznie poprzedzającej go krawędzi ( Pred[k],k ), a ta - z dokładnością do orien- tacji odpowiedniemu łukowi w grafie G. Zatem każdy węzeł k nie będący korzeniem reprezentuje dokładnie jeden łuk bazowy (zmienną bazową). Dla wygody przyjmiemy następującą konwencję: jeśli a oznacza krawędź (lub odpowiadający jej łuk w grafie G) należącą do drzewa, to symbolem [a]

oznaczać będziemy węzeł, który ją reprezentuje. Odwrotnie, jeśli w oznacza węzeł, to symbol {w } oznacza łuk bazowy reprezentowany przez ten węzeł.

Aby ułatwić odczytanie orientacji łuków grafu G odpowiadających lukom drzewa T warto posłużyć się tablicą Conf Znaczenie tej tablicy jest następu- jące: jeśli krawędź (Pred[k\, k) ma orientację zgodną z odpowiadającym jej lukiem w grafie G,to Conf[k\ = 1, w przeciwnym wypadku Conf[k\ = — 1.

Dla prawidłowego działania metody sympleks opartej na prezentowa- nej tu strukturze danych niezbędne jest również prowadzenie słownika Ro, którego kolejne składowe przypisują wierzchołkom drzewa numery pozycji w bazie B , na których znajdują się reprezentowane przez te wierzchołki kolumny (łuki). Jeśli zatem k jest zwyczajnym wierzchołkiem drzewa (nie korzeniem), to Ro[k] oznacza pozycję pola NrB zawierającą numer kolumny (zmiennej bazowej) reprezentowanej przez węzeł k

Dla pełnej charakterystyki drzewa niezbędne jest jeszcze wskazanie, któ- re węzły są korzeniami pomocniczymi. Służy do tego tablica Er. Er[w] = 0 jeśli w nie jest korzeniem. Jeśli w jest korzeniem pomocniczym wtedy Er[w]

jest numerem kolumny jądra D odpowiadającej węzłowi w (a więc Er[w] ^

°). Bardzo ważnymi z punktu widzenia rozwiązywania układów (6.1) są ta-

blice Thread i Rethread. Tablice te definiują linię obiegową w drzewie da-

Rozwiązywanie zadań programowania liniowego 23 jąc możliwość przeglądania jego węzłów w określonym porządku liniowym, oraz w porządku odwrotnym. Dla ich określenia założymy, że następniki każdego wierzchołka są uporządkowane liniowo, od lewego do prawego. W y- korzystamy prefiksową metodę przechodzenia drzewa (metoda preorder [1]).

Przechodzenie drzewa tą metodą zdefiniowane jest indukcyjnie:

1. Odwiedź korzeń r drzewa.

2. Przejdź metodą preorder poddrzewa, których korzeniami są następ- niki węzła r , według obowiązującego dla nich porządku.

Jeśli za pomocą tej metody wierzchołki drzewa odwiedzone zostały w kolejności no, tą, . . . , um, to przyjmujemy:

Thread[vi] =

Rethreacl[vm+i-i\ - um_ t-

<def

dla i = 0 , 1 , . . . , m, gdzie vq = um+i = 0.

Tablica Depth przyporządkowuje każdemu węzłowi k odległość tego wę- zła od korzenia 0 mierzoną liczbą łuków przebytych na drodze od korzenia głównego do węzła k.

Z każdym węzłem k drzewa wiążemy poddrzewo T(k), którego korzeniem jest węzeł k i które obejmuje wszystkie węzły łączące się z korzeniem 0 poprzez drogi zawierające węzeł k.

Ostatnia tablica dotyczy tych poddrzew, a mianowicie

Last[k] = numer ostatniego na linii Thread węzła w poddrzewie T(fc), k = 0 , 1 , . .. , m. Jeśli Last[k] = /, to / £ T(k ), natomiast Thread[l] £ T(k).

7. R ozw iązyw an ie układów B u a: = / i

B

= e. Przechodząc drzewo wzdłuż linii Thread określamy dwa ciągi: ciąg kolejno przebytych wierzchołków W2 , . . . , wm odpowiadających równaniom układu B n i = / , oraz ciąg {uą}, { ^ 2}, • • •, {^m} reprezentowanych przez nie łuków bazo- wych odpowiadających zmiennym bazowym ...,

Jeśli zmienimy kolejność równań i niewiadomych układu dostosowując ją do porządków wyznaczonych przez powyższe ciągi, to otrzymamy równo-

ważny układ trójkątny górny. Układ ten można łatwo rozwiązać wykorzy- stując przedstawioną uprzednio strukturę drzewa. Przedstawimy teraz trzy algorytmy rozwiązywania układu B u a; = / . Pierwszy z nich dla przypadku gdy / jest dowolnym wektorem, drugi - gdy / ma dokładnie jeden element różny od zera (+1 lub —1), a trzeci - gdy / ma dwa elementy różne od zera (+1 i -1 ).

A l . (0) Inicjalizacja X{py = f p dla węzłów p nie będących korzeniami;

k := Last[ 0];

(1) Dopóki k nie jest korzeniem głównym ( k ± 0) realizuj punkty

(a), (b) i (c):

(a) Jeśli q := Pred[k] nie jest korzeniem, to

*^{g} •— ^{g} "ł" 5 (b) .T{fc} := Conf[k] * x {k} ; (c) k := Rethread[k ) ;

A 2. Załóżmy, że f k = u (= +1 lub -1 ), oraz fi = 0 dla i ^ k (i = 1, . . . , m). Algorytm rozwiązania układu B

.

= / jest w tym przy- padku następujący

(0) := 0 dla i = 1 , . . . , m;

(1) X{ky := Conf[k] * u>;

q := Pred[k ];

(2) Jeśli q = 0, to STOP. Wszystkie niezerowe współrzędne x otrzy- mały właściwe wartości. Jeśli q / 0, to k := q i ponownie reali- zujemy kroki 1 i 2.

A 3 . Niech teraz f k = u, fi = —u oraz fi = 0 dla i ^ k,l. Oznaczmy symbolami Sk i Si drogi łączące korzeń 0 z węzłami k i / odpowiednio, oraz h niech oznacza węzeł, w którym te drogi łączą się. Algorytm jest teraz następujący:

(0) := 0 dla i — 1 , . . . , m ;

(1) Na każdej z dróg Sk i Si poniżej węzła h (dla węzłów k , dla których Depth[k] > Depth[h ]) stosujemy kolejno kroki 1 i 2 al- gorytmu A2 z wartościami

i u := —u odpowiednio.(Pozostałe współrzędne x mają wartość 0.)

Dla wyznaczenia węzła h , w którym łączą się drogi Sk i Si można wy- korzystać tablicę Depth , jak to wynika z sugestii zawartej w punkcie 1 algorytmu A3.

Koszt algorytmu A l jest proporcjonalny do liczby węzłów w sieci. Koszt algorytmu A2 (jeśli wyłączymy koszt wyzerowania wektora x) jest propor- cjonalny do długości drogi Sk w drzewie. W wielu przypadkach większość węzłów nie jest w ogóle odwiedzana w czasie działania algorytmu A2. (To samo, choć w nieco mniejszym stopniu, odnosi się do algorytmu A3.) W al- gorytmie A3 liczba działań jest proporcjonalna do sumy długości obu dróg Sk i Si pomniejszonej o długość części wspólnej tych dróg. (Również z wy- łączeniem kosztu zerowania wektora x.) Wydaje się, że osiągnięty tu został kres liczby niezbędnych dla rozwiązania układu działań i że dokonanie tego mniejszym kosztem jest niemożliwe.

Układ w B n = e jest układem o macierzy transponowanej B ^ . Jeśli

zatem uporządkujemy jego równania w kolejności łuków na linii Thread , a

niewiadome w kolejności występujących na tej linii węzłów, to otrzymamy

równoważny układ trójkątny dolny. Ten układ także da się łatwo rozwiązać

Rozwiązywanie zadań programowania liniowego 25 przy pomocy struktury drzewa. Rozważymy przy tym dwa przypadki: gdy wektor e jest dowolny i gdy ma jeden element różny od zera.

B l. Załóżmy, że e jest dowolnym wektorem.

(0) Inicjalizacja

Wk := 0 dla węzłów będących korzeniami;

i := Thread[ 0];

l := Xas<[0]];

(1) q := Pred[i]‘,

Wi := wq + * Conf[i]\

(2) Jeśli i = /, STOP; wszystkie niezerowe współrzędne w są już określone. Jeśli i / l to i := Thread[i] i powtarzamy powyższe czynności od punktu 1.

B 2. Niech e ^ } będzie jedyną współrzędną wektora e różną od zera. Zgod- nie z przyjętą umową luk {k } reprezentowany jest w drzewie przez kra- wędź (Pred[k], k) oraz węzeł k. Algorytm rozwiązania układu tuBn = e jest w tym przypadku następujący:

(0) Inicjalizacja

Wi := 0 dla każdego węzła sieci;

oj := e{k} * Conf[k}\

i := k ;

1

Wi

u;

(2) Jeśli i = Last[k ], STOP; wszystkie niezerowe współrzędne w są już określone. Jeśli i / Last[k], to i := Thread[i] i powtarzamy powyższe czynności od punktu 1.

Jeśli nie liczyć początkowego zerowania wektora w , to koszt algorytmu B2 jest proporcjonalny do liczby węzłów w poddrzewie o korzeniu w wierz- chołku k , gdzie ^ 0. Koszt algorytmu B l jest liniowo zależny od m.

Przedstawione tu algorytmy mają następującą bardzo dobrą własność:

wartość nadana dowolnej niewiadomej w dowolnym kroku algorytmu jest dokładną wartością tej niewiadomej i nie ulega zmianie w żadnym z następ- nych kroków algorytmu.

8. O dn ow ien ie drzew a przy zm ianie bazy. Zajmiemy się teraz mo-

dyfikacją drzewa T polegającą na wstawieniu do niego nowej krawędzi z

jednoczesnym usunięciem innej krawędzi leżącej na powstałym w ten spo-

sób cyklu. Krawędź wstawiana odpowiada kolumnie wprowadzanej do bazy

w jednej iteracji metody sympleksowej, a krawędź usuwana kolumnie wy-

rzucanej z bazy. Po tej wymianie zmodyfikowane drzewo T ma być nadal

drzewem rozpinającym, a więc zarówno krawędź wstawiana jak i usuwana

muszą leżeć na tym samym cyklu.

Niech eu oznacza pewną krawędź drzewa T , a ed inną. krawędź grafu G°, nie należącą do T .

Jeśli do drzewa T dołączymy krawędź ed, powstanie w nim cykl C (a mówiąc ściślej - w grafie G0), składający się z dwóch gałęzi wychodzących z jednego wierzchołka, połączonych krawędzią ed. Oznaczmy ten wspólny wierzchołek symbolem v (ma on najmniejszą głębokość ze wszystkich wierz- chołków cyklu C ). Załóżmy, że krawędź eu leży na cyklu C . Usunięcie jej spowoduje, że znów powstanie drzewo. Nazwijmy ciąg krawędzi i wierzchoł- ków cyklu C nie przechodzący przez Vh i zawarty między krawędziami eu oraz ed - ścieżką krytyczną.

Zamieszczony na rys. 2 szkic cyklu ilustruje wprowadzone pojęcia. Dołą- czenie krawędzi ed oraz usunięcie krawędzi eu powoduje, że należy odwrócić kolejność występowania wierzchołków ze ścieżki krytycznej (w sensie po- rządku Pred ), a to powoduje konieczność zmian innych funkcji opisujących drzewo.

Mając zatem dane drzewo T oraz krawędzie ed oraz eu G T takie, że eu leży na cyklu powstałym w drzewie T po dołączeniu krawędzi ed, należy przekształcić drzewo T w drzewo T' tak, aby:

1. ed G T1 oraz eu ^ T7;

2. funkcja Pred w węzłach poza ścieżką krytyczną nie uległa zmianie, natomiast na ścieżce krytycznej, w drzewie T' Vi = Pred(vi+ 1 ), i = 1 , . . . , Ar — 1;

3. funkcje pozostałe były zgodne z funkcją Pred w nowym drzewie T'.

Dla tego przekształcenia można wykorzystać algorytm, którego idea jest następująca: w kolejnych krokach odłączamy poddrzewa o korzeniach rą, V2,...,V k i podłączamy je jedno pod drugim, poniżej krawędzi ed. Podłą- czamy je zawsze jako skrajne lewe poddrzewa, co zwiększa szansę, że ojciec podłączanego poddrzewa będzie miał zachowaną wartość funkcji Last.

Algorytm składa się z trzech etapów:

I. Przełączanie kolejnych poddrzew o korzeniach v\, vo , . .. , Ufc- W trak-

cie tego przełączania ustalone zostają wartości funkcji Pred, Thread,

Rozwiązywanie zadań programowania linioiuego 27

Rethread i Depth , zaś wartości funkcji Last są modyfikowane; osta- teczne ich ustalenie następuje w II i III etapie.

II. Ostateczne ustalenie wartości funkcji Last dla wierzchołków . .. , tą.

III. Ostateczne ustalenie wartości funkcji Last dla wierzchołków leżących powyżej vq i powyżej Vk+i (w relacji drzewa T). (Szczegółowy opis algorytmu: Jabłonowski, Zorychta [7].)

9. O p era cje z ją d re m . W obu przypadkach, zarówno przy wyznacza- niu cen dualnych jak i przy obliczaniu rozkładu kolumny wchodzącej do bazy niezbędna jest znajomość macierzy V -1 lub pewnej jej reprezentacji umożliwiającej obliczanie iloczynów e V -1 lub V -1 / . Macierz tę, lub jej re- prezentację, nazywać będziemy jądrem. Zwróćmy uwagę, że jądro stanowi

” prawy dolny” (umownie, z dokładnością do kolejności wierszy i kolumn) narożnik macierzy B _1. Jeśli zechcemy korzystać z rozszerzonej macierzy V dla automatycznego wyznaczania W2 , to V -1 będzie prawym dolnym narożnikiem odwrotności rozszerzonej bazy

B = B 0

c b

1

Jednym ze sposobów prowadzenia i odnawiania jądra jest eliminacja Gaussa-Jordana na B -1 obcięta do pola V -1 . Nie jest to sposób nume- rycznie stabilny, ale za to prosty w realizacji. Ograniczenie się do zadań, w których wymiar jądra jest stosunkowo niewielki, oraz możliwość prowadze- nia obliczeń w podwójnej precyzji daje nadzieję na uzyskanie użytecznego algorytmu dla dość obszernej klasy zadań praktycznych.

Dla wygody implementacji podamy górne ograniczenie wymiaru jądra w zależności od liczby wierszy i kolumn niesieciowych w macierzy A .

L

3. Jeśli macierz A ma p niesieciowych wierszy i q niesieciowych kolumn , to stopień jądra nie przekracza p + q.

Przedstawimy teraz pewien projekt przechowywania i przekształcania jądra oraz wykonywania lewo- i prawostronnych mnożeń z jego udziałem.

Jądro zmienia się z iteracji na iterację przy czym charakterystyczne cechy tych zmian są następujące:

1. wymiar jądra może wzrosnąć lub zmaleć o 1, lub też może nie ulec zmianie;

2. zestaw wierszy i kolumn B reprezentowanych w V -1 może się zmie- nić nawet wtedy, gdy nie ulega zmianie wymiar tej macierzy;

3. zmienia się na ogół uporządkowanie wierszy i kolumn V -1 i w związku z tym jest ono przeważnie inne niż uporządkowanie standardowe.

Jądro przechowywane będzie w tablicy D o stałym wymiarze h, ograni-

czonym z góry przez p-f- q na podstawie lematu 3. W praktyce przyjmujemy

h = p + q + 1, gdyż nie można z góry przewidzieć jaki w toku procesu iteracyjnego zdarzy się maksymalny wymiar V -1 .

Zasady obsługi i wykorzystania tablicy D są następujące:

1. Wiersze i kolumny D numerowane są od 1 do h.

2. Wiersze (kolumny) bieżącej macierzy V -1 przechowywane są w od- powiednich wierszach (kolumnach) tablicy D , na ogół ze zmianą kolejności i bez zachowania ciągłości numeracji (pewne wiersze i kolumny D , nieko- niecznie końcowe, mogą być niewykorzystane.

3. Do obsługi wierszy D służy pole indeksów Rox o podanym niżej znaczeniu

rii E l,k - numer pozycji bazowej,

której odpowiada i-ty wiersz D;

— 1 - jeśli wiersze i , . . . , h są niewykorzystane;

0 - dla pozostałych i.

Tak więc jeśli Rox[i\ > 0, to NrB[rii\ jest numerem kolumny (zmiennej) ba- zowej, której odpowiada wiersz macierzy V -1 przechowywany w i-tym wier- szu tablicy D . Jeśli Rox[i] < 0, to i-ty wiersz D nie zawiera żadnego wiersza aktualnej V -1 (jest niewykorzystany, wolny). Rox[i] = — 1 jest znacznikiem końca wykorzystania wierszy D (wszystkie wiersze o numerach i , ... ,h są wolne, niezależnie od wierszy wolnych, które mogą wystąpić wewnątrz bloku zawierającego V -1 ).

4. Do obsługi kolumn D służy pole indeksów Eq o następującym zna- czeniu:

ej E l,fc - numer wiersza macierzy A ,

któremu odpowiada ji-ta kolumna D;

— 1 - jeśli kolumny j , . .. , h są niewykorzystane;

0 - dla pozostałych j.

Eq[j] > 0 wskazują bądź na równania niesieciowe bądź na równania sieciowe odpowiadające sieci transportowej, które w bazowym drzewie rozpinającym mają status korzeni pomocniczych. Jeśli Eq[j] < 0, to j - ta kolumna D jest wolna (niewyl^rzystana). Eq[j] = — 1 jest znacznikiem końca wykorzystania kolumn D (wszystkie kolumny D o numerach j , . .. , h są wolne, niezależnie od kolumn wolnych, które mogą występować wewnątrz bloku zawierającego V - 1).

5. Wiersze i kolumny scharakteryzowane są przez odpowiadające im parametry freerow i freecol , których znaczenie jest następujące:

freerow - numer pierwszego wolnego wiersza w D, freecol - numer pierwszej wolnej kolumny w D.

Parametry te mogą wskazywać bądź pozycje wewnątrz bloku zawierającego V -1 , bądź też pierwsze pozycje po końcu bloku, tuż przed znacznikiem — 1.

Eq[j] = <

Rox[i\ = <

Rozwiązywanie zadań programowania liniowego 29 Na rys. 2 zilustrowane jest rozmieszczenie w tablicy D pewnej macierzy V -1 stopnia 4.

6. Jeśli z D usuwamy wiersz i oraz kolumnę j , to na odpowiadające im pozycje w Rox i Eq wstawiamy 0. Zmieniamy przy tym wartości freerow i freecol , jeśli zwolniony wiersz lub koluma okazują się być wcześniejsze niż na to wskazują wartości tych parametrów. Jeśli po wykonaniu tych czynności znacznik końca (wierszy lub kolumn) poprzedzany jest w odpowiednim polu ( Rox lub Eq ) przez 0, to przesuwamy go o jedną pozycję wstecz.

7. Jeśli do D włączone mają być wiersz i kolumna, to wstawiamy je na pozycje wskazane przez freerow i freecol , nadając jednocześnie tym pa- rametrom nowe, zaktualizowane wartości, i przesuwając (jeśli to konieczne) znaczniki końca o 1 pozycję w przód.

8. Niezależnie od zmian wymiaru V -1 , wartości Rox[i] i Eq[j] muszą być aktualizowane zawsze gdy dochodzi do zmian składu i kolejności wierszy w D.

3 1 2 3 4 5 6 7 8

Eq 19 20 6 7 0 - 1

i Rox

1 0

2 7

3 56

4 0

5 40

6 0

7 9

8 - 1

111 112 113 114

121 122 123 124

131 132 133 134

141 142 143 144

freerow= 1, freecol— 5

Rys. 3. Ilustracja struktury jądra

Na tym tle przedstawimy kilka uwag charakteryzujących mnożenia z udziałem jądra oraz operację odnawiania jądra przy zmianie bazy. Mno- żenia te są dwóch typów: y = Da oraz ta = eD, przy czym każde z nich obcięte jest jedynie do części zawierającej V -1 (dodatkowy wiersz w D za- wiera część w2 wektora dualnego w, por. (5.6)).

Mnożenie y = Da. Mnożenia tego typu są zawsze realizowane przy zało- żeniu, że a jest częścią wektora a uporządkowanego standardowo [1,..., k], a y ma być wstawiony na odpowiednie pozycje innego wektora y o porządku standardowym [B I , . . . , Bk].

Przeszukujemy kolejno elementy pola Rox mnożąc wiersz Dj przez a jedynie wtedy, gdy Rox[i] > 0. Iloczyn wstawiamy na pozycję Rox[i] pola

y, tzn.

|/[i2oa:[t]] = D ta

Sygnałem zakończenia obliczeń wszystkich iloczynów wierszy D przez a jest napotkanie znacznika kresu wierszy V -1 , Rox[i ] = — 1.

Analogicznie postępujemy przy obliczaniu każdego iloczynu D ;«: prze- glądamy kolejno elementy pola Eq wykonując mnożenia d[i,j] * d[Eq[j]\ je- dynie dla Eq[j] > 0. Sygnałem kompletosci obliczanego sumoiloczynu jest napotkanie znacznika kresu kolumn D, Eq[j] = — 1.

Mnożenie w = eD. Symetryczną sytuację mamy przy obliczaniu iloczynu w = eD. Wektor e jest częścią wektora pełnego e uporządkowanego stan- dardowo [B I , . . . , Bk], a wynik w - częścią wektora dualnego w o porządku standardowym [1,..., k].

Dla obliczenia iloczynu w = eD przeszukujemy kolejno elementy pola Eq obliczając sumoiloczyny eD-7 (gdzie D-7 - j - ta kolumna D ) jedynie wtedy, gdy Eq[j] > 0. Wartość obliczonego sumoiloczynu wstawiamy na pozycję Eq[j] pola w, tzn.

w[Eq[j\] = eD-7.

Napotkanie znacznika kresu kolumn V -1 sygnalizuje zakończenie obliczania sumoiloczynów.

Przy obliczaniu każego sumoiloczynu eD-7 przeglądamy kolejno elementy pola Rox wykonując mnożenia e[Rox[i]] *d[i,j] jedynie wtedy, gdy Rox[i] >

0. Sygnałem kompletności obliczenia sumoiloczynu jest napotkanie znacz- nika kresu wierszy D .

Odnawianie jądra. Przyjmujemy, że odnawianie jądra dokonywane bę- dzie przy użyciu metody eliminacji Gaussa-Jordana. Zmiana jądra podzie- lona jest na dwa etapy.

I. Etap pierwszy polega na usunięciu z D wierszy i kolumn, które w następnej iteracji nie będą wchodziły w jego skład i wprowadzeniu nowych wierszy i kolumn, które będą do D należały począwszy od następnej ite- racji. Zmiany te ograniczają się do 0-2 wierszy i kolumn i zależą od typu przekształcenia (por. rozdz. 10).

Usuwanie wierszy i kolumn sprowadza się do podstawienia Rox[i\ = 0 i Eq[j] = 0 na odpowiednich pozycjach.

Nowy wiersz B -1 (oznaczmy go ft) wprowadzany jest zawsze na pozy- cję freerow przez umieszczenie jego współczynników zgodnie z kolejnością wyznaczoną przez wskaźniki Eq[j] > 0, tzn. d[freerow, j] = ft[Eq[j]] oraz podstawienie Rox[freerow] = i.

Analogicznie, nowakolumnaz B -1 (oznaczmy ją przez f t ) wprowadzana jest do D na pozycję freecol zgodnie z formułą d[i, freecol] = ft[Rox[i]]

dla Rox[i] > 0. Towarzyszy temu podstawienie Eq[freecol] = j.

Każdorazowo po dołączeniu lub usunięciu wiersza lub kolumny do lub z

D modyfikowana jest odpwiednia wartość freerow lub freecol oraz (ewn-

tualnie) znacznik końca Rox[i\ = -1 lub Eq[j] = — 1.

Rozwiązywanie zadań programowania linioiuego 31 II. Etap drugi to proces eliminacji na zmodyfikowanej w etapie pierw- szym tablicy D . Do zrealizowania tej operacji niezbędne są następujące dane:

- rozkład ćf kolumny a wprowadzanej do bazy (por.(5.7)), - numer r wiersza głównego,

- część wiersza f3r macierzy B -1 przypadająca na jądro.

Potrzebna część wiersza j3r może być w pewnych sytuacjach dostępna bezpośrednio jako jeden z wierszy D lub może być wyznaczona z użyciem wzoru (5.9).

Przyjmijmy oznaczenie j3'r = o " 1/?,.. Wiersze tablicy D, dla których Rox[i] > 0 oraz Rox[i] ^ r przekształcane są następująco:

D ' = Dj - a[Rox[i]]f3'r.

Jeżeli wiersz główny ma należeć do D, to został on już włączony do D w etapie I i teraz należy go jedynie zastąpić przez fi'r podstawiając jednocześnie Rox[g\ = r, gdzie g jest numerem wiersza w D , w który wstawiony został wiersz f3'r.

10. Typy przekształceń bazy. Jeśli w danej iteracji metody sym- pleks wybrana już została zmienna x^s wchodząca do bazy i zmienna a*Br opuszczająca bazę, wtedy należy zmienić reprezentację bieżącej bazy, a w szczególności odnowić bazowe drzewo rozpinające oraz jądro. Zmiana drzewa polega na odnowienu wartości opisujących go funkcji, a zmiana jądra na ad- aptacji tablicy D i związanych z nią parametrów pomocniczych.

Wymiana luku w drzewie odbywa się na jednej z następujących zasad:

1) Łuk a £ B2 można przesunąć do B I bez usunięcia innego luku z B I jeśli pętla zamknięta w drzewie przez przesuwany luk a zawiera luk pomocniczy b. Należy wtedy

- wprowadzić do drzewa luk a,

- zdjąć z węzła [b] reprezentującego luk b status korzenia pomocniczego i usunąć luk pomocniczy b ,

- usunąć z jądra wiersz odpowiadający łukowi b i kolumnę odpowiadającą węzłowi [b].

2) Łuk a £ B I można przesunąć z B I do B 2 bez wprowadzania w jego miejsce innego luku, w następujący sposób

- usunąć z drzewa luk o,

- nadać węzłowi [«] reprezenującemu dotąd luk a status korzenia pomoc- niczego i podwiązać go do korzenia głównego lukiem pomocniczym,

- rozszerzyć jądro o wiersz macierzy B -1 odpowiadający lukowi a i ko- lumnę odpowiadającą węzłowi \a].

3) Łuki a £ B I i b £ B2 można zamienić miejscami na dwa sposoby:

1. jeśli oba luki a i b należą do tej samej pętli w drzewie, to można dokonać

prostej wymiany wprowadzając do drzewa luk b i usuwając a ; 2. jeśli te

luki nie należą do wspólnej pętli to można próbować stosować operacje 1) i 2) w odpowiedniej kolejności (zamiana będzie możliwa przy spełnieniu wymaganych przez te operacje warunków).

4) Dowolny luk pomocniczy a może być zastąpiony przez inny luk po- mocniczy w następujący sposób:

- w poddrzewie korzenia pomocniczego [a] wybieramy dowolny węzeł w nadając mu status korzenia pomocniczego i podwiązując do korzenia głów- nego lukiem pomocniczym,

- rozcinamy powstałą pętlę zmieniając status węzła [a] na zwykły i usu- wając z niej łuk pomocniczy a,

- w miejsce kolumny odpowiadającej węzłowi [a] wprowadzamy do D kolumnę odpowiadającą nowemu korzeniowi pomocniczemu w.

Ze zmianami drzewa związane są ściśle zmiany jądra D. 0 zmianach współczynników D , które są konsekwencją stosowania metody eliminacji była już mowa w rozdz. 7. Oprócz tych zmian istotne są zmiany strukturalne, polegające na zmianie zestawu wierszy i kolumn B-1 wchodzących w jądro D . Klasyfikacja tych zmian jest następująca:

1. bez zmian;

2. dodanie wiersza w miejsce innego, usuniętego wiersza;

3. dodanie kolumny w miejsce innej, usuniętej kolumny;

4. usunięcie wiersza i kolumny;

5. dodanie wiersza i kolumny.

Szczegółowa analiza zasad wymiany w drzewie i w jądrze, jak również mechanizmu działania opisanego w rozdz. 3 algorytmu sympleks pozwala wyliczyć siedem typów przekształceń bazy. Dla wygody prezentacji tych typów przyjmiemy następujące oznaczenia:

NET - zbiór zmiennych lub kolumn sieciowych;

LOOP(a;t) - pętla w drzewie powstała w wyniku dostawienia dodat- kowego łuku a odpowiadającego zmiennej Xi (równoważne oznaczenie - LOOP(a));

EA - łuk pomocniczy.

Typ 1. 1. xbt G XB2

2. a?Ns £ NET lub jeśli a*NS G NET to EA ^ LOOP(

s) Operacje :

1. eliminacja (z wierszem głównym r należącym do D ).

Typ 2. 1. x bt G XB2 2. a?NS G NET 3. EA <E LO OP(

Ns) Operacje :

1. wstawienie do drzewa łuku x^s i usunięcie EA z pętli LOOP(:

)

jednoczesnym pozbawieniem węzła [EA] statusu korzenia pomocniczego;

2. usunięcie z D kolumny odpowiadającej [EA];

Rozwiązywanie zadań programowania liniowego 33 3. eliminacja (z wierszem głównym r należącym do D );

4. usunięcie po eliminacji wiersza głównego r.

Typ 3. 1. xBr G a?Bi

2. istnieje x t G xB2 taki, że xBr G LOOP(:r*)

3. #Ns ^ NET lub jeśli #

G NET to ^Br ^ LOOP(#

s) oraz EA ^ LOOP(

s)

Operacje:

1. usunięcie z D wiersza odpowiadającego zmiennej x t ; 2. dołączenie do D wiersza r macierzy B -1 ;

3. eliminacja (z wierszem głównym r włączonym do D );

4. usunięcie z drzewa luku a?Br i wprowadzenie luku xt.

Typ 4. 1. xBr G x b i

2. nie istnieje x t G xB2 taki, że xBr € LOOT(xt)

3. #Ns $ NET lub jeśli zNs G NET to zBr i LOOP(zNs) oraz EA ^ LOOP(£ n s)

Operacje:

1. dołączenie do D wiersza r macierzy B -1 ;

2. dołączenie do D kolumny / macierzy B -1 , gdzie / = Last[rl] i Ro[r'] = r;

3. eliminacja (z wierszem głównym r włączonym do D );

4. usunięcie z drzewa luku

f , nadanie węzłowi / statusu węzła pomoc- niczego i podwiązanie go do korzenia głównego lukiem pomocniczym.

Typ 5. 1. z Br £ x b i 2. #

G NET 3. a?Br E LOOP(:rNs)

4. jeśli EA G LOOP(

s) to nie istnieje x t G xB2 taki, że xBr G LOOF(xt)

Operacje:

1. wyznaczenie wiersza głównego eliminacji (r-ty wiersz macierzy B _1);

2. eliminacja;

3. wprowadzenie do drzewa luku £

i usunięcie luku XBr

#Br G

^Ns £ NET EA G LOOP(

)

istnieje x t G xB2 taki, że xBr G LOOP(x^) TypQ. 1.

2 ^. 3.

Operacje: 4.

1. usunięcie z D wiersza odpowiadającego zmiennej xf, 2. usunięcie z D kolumny odpowiadającej węzłowi [EA];

3. wyznaczenie wiersza głównego eliminacji (wiersz r-ty macierzy B 4. eliminacja;

5. wprowadzenie do drzewa luku

i usunięcie luku pomocniczego EA z jednoczesnym pozbawieniem węzła [EA] statusu korzenia pomocniczego;

! — 1 _);

6. wprowadzenie do drzewa łuku xt i usunięcie łuku XBr*

Typ 7. 1. a?Br G #

2. #Ns £ NET 3. EA

LOOP(

s) 4. xBr £ LOOP(

Ns)

5. nie istnieje xt G £

2 taki, że xBr G LOOP(xt) Operacje:

1. usunięcie z D kolumny odpowiadającej węzłowi [EA];

2. dołączenie do D kolumny l macierzy B -1 , gdzie / = Last[r'] i Ro[r '] = r;

3. wyznaczenie wiersza głównego eliminacji (wiersz r-ty B -1 );

4. eliminacja;

5. nadanie węzłowi l statusu korzenia pomocniczego i podwiązanie go do korzenia głównego lukiem pomocniczym z jednoczesnym usunięciem z drzewa łuku a?Br5

6. zdjęcie z węzła [EA] atrybutu korzenia pomocniczego, usunięcie łuku EA z jednoczesnym wprowadzeniem do drzewa łuku #

s.

Wyliczone typy przekształceń obejmują wszystkie możliwe przypadki ja- kie mogą się zdarzyć przy zmianie bazy. Oczywiście, w każdej iteracji al- gorytmu sympleks realizowane jest tylko jedno z powyższych przekształceń.

Dla przekształceń tych spełniony jest następujący lemat:

L

4. Jeżeli w bieżącej iteracji stopień B u jest maksymalny , to po przekształceniu bazy za pomocą jednego z powyższych przekształceń stopień nowego B u będzie nadal maksymalny.

11. Zakończenie. Aby zrealizować algorytm sympleks opisany w roz- dziale 3 wystarczy w odpowiednich jego krokach wykorzystać struktury i algorytmy przedstawione w kolejnych rozdziałach pracy (5-10). Dotyczy to kroku 1, gdzie obliczamy wektor dualny w (3.2), kroku 3 w którym obliczamy rozkład a (3.6) kolumny wchodzącej do bazy, i kroku 5, kiedy to struktura bazy (drzewo + jądro) jest odnawiana po wymianie kolumny w bazie.

Proces można zainicjować startując ze sztucznej bazy jednostkowej. W części sieciowej odpowiada to dołączeniu m sztucznych węzłów, z których każdy podłączony jest bezpośrednio do korzenia głównego. W kolejnych ite- racjach sztuczne luki są z bazy usuwane. Odpowiada to usuwaniu z drzewa sztucznych węzłów i łuków i wprowadzaniu na ich miejsce węzłów i łuków sieci transportowej. Sztuczny węzeł (wraz z odpowiadającym mu lukiem), który został z drzewa usunięty, nigdy już do drzewa nie wraca.

Rozmiar pamięci potrzebnej dla pomieszczenia struktur opisujących

część sieciową bazy zależy liniowo od liczby węzłów sieci. Również czas

wykonania operacji z udziałem drzewa (mnożenie lewe i prawe, odnawia-

nie drzewa) jest proporcjonalny do liczby węzłów. Stąd, przy założeniu że

Rozwiązywanie zadań programowania liniowego 35 wymiar jądra jest mały w stosunku do wymiaru drzewa, koszt jednej itera- cji sympleksowej jest dla dużych zadań znacznie mniejszy niż 0 (k 2). Warto podkreślić również to, że z uwagi na mniejszą, niezbędną do obliczeń pamięć możliwe jest rozwiązywanie względnie dużych zadań na mikrokomputerach typu IBM P C /X T , a więc jedynym, w miarę powszechnie dostępnym u nas sprzęcie.

Opisana w pracy metoda SON rozwiązywania zadań programowania li- niowego z wbudowaną strukturą sieciową została zastosowana w systemie DINAS (Ogryczak, Studziński, Zorychta [8,9]) służącym do analizy wielo- kryterialnych zadań transportowo-lokalizacyjnych. Dotychczasowe doświad- czenia z tym systemem potwierdzają skuteczność metody SON w rozwiązy- waniu względnie dużych zadań o opisanej strukturze na sprzęcie typu mikro.

Literatura

[1] A. V. A h o, J. E. H o p c r o ft, J. D. Ulm an, The design and analysis of computer algorithms, Addison-Wesley 1974.

[2] G. B. D a n tz ig , Linear programming and extensions, Princeton University Press 1963.

[3] F. G lo v e r , D. K lin g m a n , The simplex SON method for LPembedded network problems, Mathematical Programming Study 15, 1981.

[4] F. G lo v e r , D. K lin g m a n , Basis exchange characterizations for the simplex SON algorithms for LPembedded networks, Mathematical Programming Study 24, 1985.

[5] M. D. G r ig o r ia d is , An efficient implementation of the network simplex method, Mathematical Programming Study 26, 1986.

[6] T. C. Hu, Integer programming and network flows, Addison-Wesley 1970.

[7] J. J a b ło n o w s k i, K. Z o r y c h ta , Pierwotna metoda sympleksowa dla przepływów sieciowych, Matematyka Stosowana XXXII, 1989.

[8] W. O g r y c z a k , K. S tu d z iń sk i, K. Z o ry ch ta , A solver for the multi-objective transshipment problem with facility location, European Journal of Operational Re- search 43, 1989. (również: Sprawozdania Instytutu Informatyki Uniwersytetu War- szawskiego, Nr 169, Wyd. Uniw. Warsz., 1988)

[9] W. O g ry cz a k , K. S tu d z iń sk i, K. Z o ry ch ta , General concepts of the Dynamie Interactive Network Analysis System (DINAS), Archiwum Automatyki Telemecha- niki, t. XXXII, z. 4, 1987.

[10] K. Z o r y c h t a , W. O g ry cz a k , Programowanie liniowe i całkowitoliczbowe, Wydaw- nictwa, Naukowo-Techniczne, Warszawa 1982.

Summary

On solving linear programming problems with embedded network structure

A special partitioning algorithm for solving linear programming problems with embed-

ded network structure is presented. As an example of such a problem the minimum-cost

network flow problem under additional linear constraints can be considered. This algori-

thm is a primal simplex basis partitioning method that uses special updating and labeling

procedures to accelerate computations involving the network linear programming inter-

face. These procedures are discribed in detail to develop an efficient implementation of

the method.