papír pro razítko ústavu K611

papír pro razítko ústavu K611

1

České Vysoké Učení Technické v Praze Fakulta Dopravní

BAKALÁŘSKÁ PRÁCE

Miroslav Vaniš Matematické modelování vybraných problémů v dopravě v jazyce Java Ústav Aplikované Matematiky

Vedoucí bakalářské práce: Studijní program: Studijní obor:

Mgr. Martin Scholtz, Ph.D. Technika a technologie v dopravě a spojích Automatizace a informatika

Praha 2013

Papír se zadáním Bc. práce

Poděkování Na tomto místě bych rád poděkoval všem, kteří mi poskytli podklady a cenné rady pro vypracování této práce. Zvláště pak děkuji Mgr. Martinu Scholtzovi, PhD. za odborné vedení a konzultování bakalářské práce i za rady, které mi poskytoval během celého mého studia.

Prohlášení Prohlašuji, že jsem tuto bakalářskou práci vypracoval(a) samostatně a výhradně s použitím citovaných pramenů, literatury a dalších odborných zdrojů a za pomoci svého školitele. Prohlašuji, že jsem predloženou práci vypracoval samostatne a že jsem uvedl veškeré použité informacní zdroje v souladu s Metodickým pokynem o etické príprave vysokoškolských záverecných prací. Nemám závažný duvod proti užívání tohoto školního díla ve smyslu 60 Zákona c.121/2000 Sb., o právu autorském, o právech souvisejících s právem autorským a o zmene nekterých zákonu (autorský zákon).

V Praze dne . . . . . . . . . . . .

Podpis autora

Abstract In this bachelor thesis we have studied selected mathematical models developed to mimic the behavior of the real traffic flow. In the theoretical part we introduced some basic ideas related to the models of the traffic flow. We introduce both macroscopic and microscopic models, then we restrict our attention to the latter and briefly review existing microscopic models. Next we introduce the elementary theory of dynamical systems and discuss the stability analysis by linearization. Finally we present numerical methods for solving ordinary differential equations (Euler method, Runge-Kutta methods) and present the so-called Intelligent Driver Model (IDM) in detail. In the practical part of the thesis we implement the IDM model in computer language Java. Our application, being the most important result of the thesis, solves dynamical equations of IDM model and represents the results both in the graphic form and in the form of data files. Moreover, we have extended basic IDM model by overtaking model known as MOBIL. At the end of the thesis, after explaining main features and algorithms used in the application, we present full source code and a number of results obtained.

Abstrakt V této bakalářské práci studujeme vybrané matematické modely, které se snaží napodobit chování dopravního toku. V teoretické části představíme základní principy spojené s těmito modely. Uvedeme jak makroskopické tak i mikroskopické modely. Poté se omezíme pouze na mikroskopické modely a krátce je popíšeme. Dále se seznámíme se základní teorií dynamických systémů a jejich analýzou stability po linearizaci. Nakonec ukážeme numerické metody pro řešení obyčejných diferenciálních rovnic (Eulerova metoda, metody Runge-Kutta) a podrobně vysvětlíme, jak funguje Intelligent Driver Model (IDM). V praktické části této práce implementujeme IDM model do programovacího jazyka Java. Tato aplikace, která je hlavním výsledkem této práce, řeší rovnice IDM modelu a zobrazuje výsledky jak v grafické podobě, tak i ve formě datových souborů. Navíc rozšíříme IDM model předjížděcím modelem MOBIL. Na konci práce, po vysvětlení základních vlastností programu, je uveden celý zdrojový kód aplikace a dosažené výsledky.

6

Obsah Úvod Cíle práce . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3 4

I

5

Teoretický úvod

1 Makroskopické modely 1.1 Rovnice kontinuity . . . . . 1.2 Pohybové rovnice kontinua 1.3 Doprava jako kontinuum . . 1.4 Fundamentální diagram . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

6 6 7 8 9

2 Mikroskopické modely

11

3 Dynamické systémy 3.1 Planární dynamický systém . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Klasifikace kritických bodů . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14 14 16

4 Numerické metody 4.1 Dynamický systém . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Eulerova metoda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Metoda Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19 19 20 21

5 Intelligent Driver Model 5.1 Matematický popis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Model MOBIL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.1 Popis modelu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23 24 24 25

II

28

Vypracování

6 Popis programu 6.1 Převod kartézských souřadnic 6.2 Balík Utilities . . . . . . . . . 6.2.1 Třída Conversions . . 6.2.2 Třída Vectors . . . . . 6.3 Balík mathcomps . . . . . . . 6.3.1 Třída JDrawPanel . . 6.4 Třída ControlPanel . . . . . . 6.5 Třída CarList . . . . . . . . . 6.6 Třída Highway . . . . . . . . 6.7 Třída FileCreator . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

1

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

29 29 30 30 31 32 32 32 33 34 35

7 Zdrojový kód 7.1 Hlavní program . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Knihovna Utilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Knihovna mathcomps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37 37 51 52

8 Simulace 8.1 Stacionární stav . . . . . . . . . . . 8.2 Náhodné rychlosti . . . . . . . . . 8.3 Změna parametrů prvního vozidla 8.4 Uzavření pruhu komunikace . . . .

55 55 57 58 62

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

Závěr

65

Literatura

67

Seznam obrázků

68

Seznam tabulek

69

2

Úvod Matematika bývá často nazývána „královnou vědÿ. Je zvykem rozlišovat mezi čistou matematikou a matematikou aplikovanou. Zatímco čistá matematika zkoumá abstraktní objekty a vztahy mezi nimi, aplikovaná matematika se snaží tyto teoretické poznatky uplatnit při řešení konkrétních reálných problémů. Možnosti aplikace matematiky jsou vskutku široké: fyzika a teoretická fyzika, chemie, teorie řízení, ekonomie, sociologie, stavebnictví, elektronika, inženýrské aplikace, dopravní problémy a mnoho jiných. Zcela speciální je pak vztah matematiky a fyziky. Pro fyziku není matematika jen nástrojem k řešení problémů, ale i jazykem, kterým jsou přírodní zákony a jevy popisovány. Není příliš překvapivé, že k analýze různých fyzikálních problémů, které se často prolínají s těmi inženýrskými, je matematika nutná a užitečná. Je zřejmé, že k navržení a postavení mostu či budovy jsou nutné přesné propočty založené na fyzikálních zákonech. Ačkoli jsou tyto zákony relativně jednoduché (Newtonovy zákony), jejich aplikace na konkrétní problém může být značně složitá a k řešení je třeba použít pokročilých matematických metod. Je třeba zdůraznit, že fyzika je ze všech vědních disciplín daleko nejvíc prostoupena matematikou a k řešení základních i aplikovaných problémů fyziky byly vyvinuty velmi efektivní a robustní matematické metody, například metody statistické fyziky, renormalizační metody, aproximační a perturbační metody, teorie integrabilních a neintegrabilních systémů, teorie stability a chaosu a mnoho jiných. Postupem času se však stále víc matematizují i jiné vědní disciplíny jako ekonomie nebo dopravní vědy. Tyto disciplíny mají společné to, že studují velmi složité systémy. To znamená, že tyto systémy sestávají z velkého počtu složek, které navíc interagují složitým způsobem. Ve fyzice se složité systémy studují v teorii tuhých látek, v molekulové fyzice, v teorii fázových přechodů a kritických jevů či v problémech mnoha těles. Proto často dochází k tomu, že problém studovaný v rámci nějaké vědní disciplíny má svůj protějšek v analogickém problému ve fyzice. Tyto analogie jsou velmi užitečné, neboť umožňují aplikovat dobře rozvinutý matematický aparát na zdánlivě zcela odlišný problém. Dopravní problémy jsou exemplárním příkladem, kde podobné analogie nacházejí využití. Dopravní systém sestává z velkého počtů účastníků, vozidel a řidičů, kteří spolu interagují složitým a ne úplně známým způsobem. Na dopravní systém lze tedy nahlížet v zásadě dvěma způsoby: jako na kontinuální prostředí, jehož dynamika je zachycena soustavou parciálních diferenciálních rovnic (obecně nelineárních), nebo jako na soustavu diskrétních elementů, přičemž dynamika každého z nich je popsána obyčejnými diferenciálními rovnicemi (obecně nelineárními). Oba způsoby popisu našly již dříve uplatnění ve fyzice v oblastech známých jako mechanika kontinua a molekulová fyzika. Diferenciální rovnice popisující dopravní nebo fyzikální problémy jsou většinou natolik složité, že je nelze analyticky řešit, a tak je nutné k jejich řešení použít numerické algoritmy. Tyto algoritmy jsou pak implementovány ve vhodném počítačovém programovacím jazyce. Inženýr, stejně tak jako fyzik, který chce počítačů využít k řešení složitých úloh se tak musí seznámit nejen se samotným problémem, ale i s celou řadou numerických metod, což je samostatná matematická oblast. Ke studiu dopravních systémů je nutná široká škála znalostí z matematiky, teorie dopravy, numerické matematiky, programování, a výhodou je i znalost fyzikálních teorií, které se svou matematickou formulací podobají problémům z dopravy. Programovat a analyzovat plnohodnotné dopravní simulace je proto možné až po dlouhém studiu všech uvedených oblastí. Cílem této bakalářské práce je učinit první krok v tomto studiu, jelikož v rámci bakalářského studia není možné zvládnout problematiku v celé její šíři. Proto cílem této práce není získání zásad-

3

ních nových výsledků, ale sestavení vlastního funkčního programu, který implementuje některý z existujících modelů.

Cíle práce Studijní cíle práce lze formulovat následovně. • Seznámení se se základními typy modelů v dopravě, tedy s modely makroskopickými a mikroskopickými. • Seznámit se s teorií dynamických systému a analýzou stability založenou na linearizaci dynamického systému. • Nastudovat základní a pokročilé numerické metody řešení soustav obyčejných diferenciálních rovnic. • Seznámit se s matematickým popisem mikroskopického modelu IDM (Intelligent Driver Model). Tato témata jsou náplní první, teoretické části práce. V druhé, praktické, části byly cíle práce následující. • Implementace modelu IDM v jazyce Java – Implementace Eulerovy metody a metody Runge-Kutta 4. řádu k řešení dynamických rovnic modelu IDM. – Využití základních grafických knihoven jazyku Java k zobrazení simulace. – Naprogramovat vlastní grafickou knihovnu, která umožní vykreslení dálnice z izometrického pohledu. – Rozšíření modelu IDM o předjíždění vozidel podle modelu MOBIL. • Analýza vlastností modelu IDM – Využít teorie dynamických systémů k analýze stacionárních stavů modelu IDM. – Implementovat metodu sečen k numerickému řešení rovnic pro stacionární stav. – Provést simulace na bázi modelů IDM a IDM+MOBIL pro stacionární stav a vybrané dopravní situace. – Makroskopický popis výsledků simulace v pojmech distribucí počtu vozidel. – Mikroskopický popis výsledků simulace. V kapitole 7 uvádíme i kompletní výpis zdrojového kódu vytvořeného programu.

4

Část I

Teoretický úvod

5

1. Makroskopické modely Makroskopické modely nahlížejí na provoz jako na spojitý dopravní proud, který je popsán spojitě se měnícími veličinami, například hustotou ρ. Makroskopické modely tedy nezkoumají skutečnou diskrétní strukturu dopravního toku, jenž je v konečném důsledku tvořen jednotlivými vozidly. Tento pohled je inspirován mechanikou kontinua, což je fyzikální disciplína, která analogicky nepřihlíží na molekulární strukturu hmoty a popisuje její rozložení na makroskopické úrovni pomocí spojitých veličin. Tyto modely jsou užitečné v situacích, kdy je hustota vozidel dostatečně velká a již poměrně malý úsek komunikace obsahuje dostatečně velký počet vozidel.[1, str.15-17] Dopravní tok se za těchto okolností chová jako spojité prostředí, jehož matematický popis je podobný popisu proudění tekutiny. V této kapitole uvedeme některé základní myšlenky, na nichž jsou makroskopické modely založeny.

1.1

Rovnice kontinuity

Jedna z nejzákladnějších rovnic mechaniky kontinua je rovnice kontinuity, která vyjadřuje zákon zachování hmotnosti. Podle tohoto zákona se hmotnost kapaliny nemůže ztrácet ani vznikat, může pouze proudit z jednoho místa do druhého. Nechť V je objem libovolně zvolené třírozmérné oblasti ohraničené uzavřenou plochou S. Celková hmotnost kapaliny uvnitř oblasti V je dána vztahem Z M = ρ dV. V

V Eulerově popisu kontinua se, na rozdíl od popisu Lagrangeova, oblast V nepohybuje spolu s kapalinou, nýbrž zůstává fixována v prostoru. Do této oblasti pak přitéká kapalina z okolí a současně vytéká z oblasti pryč. Proto se objem V nemění, ale mění se hmotnost kapaliny v tomto objemu. Změna hmotnosti kapaliny za jednotku času je dána derivací Z Z d ∂ρ dM = ρ dV = dV. (1.1) dt dt ∂t V

V

Obsahem zákona zachování hmotnosti je, že pokud se hmotnost M mění (a derivace (1.1) je nenulová), není to proto, že by nějaká hmotnost uvnitř dané oblasti vznikla či zanikla, ale proto, že z této oblasti odtekla nebo do ní přitekla. Hmotnost kapaliny, která za jednotku času proteče ohraničující plochou S je dána plošným integrálem I I Φ = ρ v · dS = j · dS, (1.2) S

S

kde j = ρv

(1.3)

je objemová hustota hybnosti kapaliny a skalární součin j · dS vyjadřuje, že do integrálu přispívá pouze složka rychlosti normálová k ploše v daném bodě. Zachování hmotnosti nyní vyžaduje, aby se integrály (1.1) a (1.2) rovnaly co do velikosti a měly různá znaménka: dM = −Φ. dt 6

Podintegrální výrazy v rovnicích (1.1) a (1.2) není možné přímo porovnat, protože první integrál je objemový a druhý plošný. Aplikace Gaussovy-Ostrogradského věty na plošný integrál pak vede ke vztahu  Z  ∂ρ + ∇ · j dV = 0. (1.4) ∂t V

Obecně platí, že je-li integrál nějaké funkce na daném intervalu roven nule, nelze z toho usoudit na nulovost integrované funkce1 . Mohlo by se například, v souladu se zákonem zachování hmotnosti, stát, že v jednom bodě prostoru vznikne hmotnost m a současně v jiném vzdáleném bodě hmotnost kapaliny klesne o tutéž hodnotu m. Pokud se nyní zvolí objem V tak, aby obsahoval oba body, pak celková změna hmotnosti je rovna nule a rovněž tak celkový tok Φ. Vztah (1.4) je tedy v tomto umělém případě splněn a zákon zachování hmotnosti platí globálně. Fyzikálně však k uvedenému procesu nemůže dojít. Proto je požadována silnější verze zákona zachování, která je lokální. To znamená, že vznikne-li v nějakém místě hmotnost dm, musela hmotnost kapaliny o dm poklesnout v bezprostředním nejbližším okolí tohoto bodu, což pouze znamená, že se hmotnost dm spojitě přemístila do blízkého bodu. Lokální zákon zachování vyžaduje, aby byla rovnost (1.4) splněna pro libovolně zvolený objem V . Dříve uvažovaný umělý případ této podmínce nevyhovuje, protože pokud je objem zvolen tak, aby obsahoval pouze jeden z bodů, bude změna hmotnosti uvnitř objemu nenulová, ale tok ohraničující plochou bude nulový. Ovšem má-li platit vztah (1.4) pro libovolně zvolený objem, musí být nule roven samotný podintegrální výraz (viz poznámka 1 pod čarou na straně 7). Rovnice kontinuity tedy zní ∂ρ + ∇ · j = 0. ∂t

1.2

(1.5)

Pohybové rovnice kontinua

Rovnice kontinuity ještě plně neurčuje proudění kapaliny (pouze v nejjednodušších případech) a je k ní potřeba dodat pohybovou rovnici. V teorii kontinua se ukazuje, že podmínkou rovnováhy kontinua je rovnice 3 X ∂σij j=1

∂xj

= 0,

kde σij je symetrický tenzor napětí. Podobně jako v mechanice soustavy hmotných bodů, i v teorii kontinua lze aplikovat d’Alembertův princip a považovat uvedenou rovnici i za pohybovou rovnici kontinua v obecném dynamickém případě. Pro izotropní prostředí je tenzor napětí diagonální a úměrný jednotkovému tenzoru, takže platí σij = −P δij , 1 Například

rovnost

Zb f (x)dx = 0. a

Je-li například funkce f antisymetrická kolem bodu (a + b)/2, je její integrál automaticky roven nule, ačkoli sama funkce f je na tomto intervalu nenulová. Ovšem, platí-li uvedená integrální rovnost na libovolném intervalu, musí být nulová i funkce f .

7

kde δij je Kroneckerovo delta (jednotková matice) a veličina P představuje tlak kapaliny. Tato volba tenzoru napětí vede na Eulerovu pohybovou rovnici ∂vi ∂vi 1 ∂P + vj =− . ∂t ∂xj ρ ∂xi V obecně anizotropním a viskózním, ale nestlačitelném prostředí je tenzor napětí dán vztahem   ∂vi ∂vj σij = −P δij + α + , ∂xj ∂xi kde α je konstanta související s viskozitou kapaliny a předpokládá se, že viskózní síly jsou úměrné gradientu rychlosti kapaliny. Je-li kapalina nestlačitelná, rovnice kontinuity se redukuje na ∇·v =0

neboli

X ∂vj = 0. ∂xj j

V tomto případě je pohybovou rovnicí kontinua Navierova-Stokesova rovnice ∂vi 1 ∂P ∂vi + vj =− + ν∆vi , ∂t ∂xj ρ ∂xi kde ν je kinematická viskozita a ∆ je Laplaceův operátor. Tato rovnice již dostatečně dobře popisuje turbulentní i laminární proudění nestlačitelné kapaliny. Obtížnost řešení NavierovýchStokesových rovnic je dobře známa a představuje již více než století výzvu pro matematiky, fyziky i experimentátory.

1.3

Doprava jako kontinuum

Z matematického hlediska vyjadřuje rovnice kontinuity nejen zákon zachování hmotnosti, ale zákon zachování jakékoliv spojité veličiny. Toho si všimli i dopravní inženýři a začali zkoumat, zda by makroskopický popis kontinua bylo možno aplikovat i na dopravní tok. V takovém případě by rovnice kontinuity vyjadřovala zákon zachování počtu vozidel. Počet vozidel samozřejmě není spojitá veličina, ale v případě vysoké hustoty dopravy je počet automobilů podstatně větší, než flukuace počtu vozidel. Nechť poloha automobilu na komunikaci je dána jedinou souřadnicí x. Zvolíme nyní délku ∆x takovou, že počet automobilů na úseku (x, x + ∆x) je dostatečně veliký; tento počet označíme n(t, x), kde jsme vyznačili, že počet vozidel závisí na poloze i na čase. Množinu automobilů, které se v tomto úseku nacházejí označíme symbolem M (t, x). Platí potom n(t, x) = |M (t, x)|, kde svislé čáry znamenají mohutnost (počet prvků) množiny. Hustotu vozidel k v bodě x a čase t definujeme vztahem k(t, x) =

n(t, x) . ∆x

(1.6)

Jinými slovy, hustota k udává počet automobilů na jednotku délky (kilometr) v bodě x. Hustota počtu vozidel k je analogon k hustotě ρ zavedené v mechanice kontinua. Takto definovanou hustotu budeme považovat za spojitou veličinu[4, str. 265]. 8

Předpokládáme, že v libovolně zvoleném úseku (x, x + ∆x) se nachází velký počet automobilů n(t, x), a každý z nich se pohybuje jinou rychlostí. Je-li však hustota toku dostatečně vysoká, takže makroskopický popis má smysl, lze očekávat, že rychlosti jednotlivých vozidel se nebudou příliš lišit. Má proto smysl definovat průměrnou rychlost toku vztahem v(t, x) =

1 X vj , n(t, x) j

kde se sčítá přes všechna vozidla z množiny M (t, x) a vj je rychlost j−tého vozidla[4, str. 274276]. Můžeme si tedy představit, že všechna vozidla na úseku (x, x + ∆x) se pohybují společnou průměrnou rychlostí v(t, x). Nakonec definujeme hustotu toku (nebo krátce tok) vztahem q(t, x) = k(t, x) v(t, x).

(1.7)

Veličina q představuje hustotu počtu automobilů, které v okamžiku t přejíždějí úsek (x, x + ∆x). Postupem podobným jako v sekci 1.1 (ale zjednodušeným na jednorozměrný případ) lze nyní odvodit rovnici kontinuity neboli rovnici zachování počtu automobilů: ∂q ∂k + = 0, ∂t ∂x

(1.8)

[1, str. 16],[4, str. 279]. Rovnice kontinuity v dopravě představuje vazbu mezi hustotou dopravního toku (počtu aut na jednotku délky) a rychlostí automobilů. Pokud si tedy dopravní inženýr zvolí určitý úsek dálnice, na kterém bude chtít měřit hlavní dopravní parametry a posléze sestavit vhodný makroskopický model, jedna z rovnic pro jeho model bude rovnice kontinuity. Jaký bude přesně zvolen model, rozhoduje právě vztah mezi hustotou a průměrnou rychlostí v oblasti. Pokud se zvolí model, kde rychlost bude záviset výlučně na hustotě v = v(ρ),

(1.9)

[4, str. 283],stačí již znát pouze rovnici kontinuity a je možno na silnici předvídat stavy budoucí, samozřejmě v závislosti na počátečních a okrajových podmínkách. Počáteční podmínkou se rozumí iniciální distribuce v čase t = 0, tedy libovolně zvolená diferencovatelná funkce k0 (x) taková, že k(0, x) = k0 (x). Vedle toho je nutné specifikovat okrajovou podmínku, tedy hustotu ρ v libovolném čase, ale pevně zvoleném bodě x0 ; bez újmy na obecnosti lze zvolit x0 = 0, takže okrajová podmínka má tvar k(0, t) = θ(t), kde θ je zase libovolně zvolená funkce času. Tyto podmínky již jednoznačně určují řešení rovnice kontinuity[4, str. 283-285].

1.4

Fundamentální diagram

V předchozí kapitole jsme vysvětlili, že k rovnici kontinuity (1.8) je nutno přidat vztah (1.9) mezi hustotou a tokem. Takový vztah se nazývá fundamentální a jeho grafická reprezentace se 9

q

v vf

kj

kj  2

k

a)

kj

k

b)

Obrázek 1.1: Fundamentální diagramy pro Greenshieldsův model. a) Závislost v = v(k). b) Závislost q = q(k). nazývá fundamentální diagram. Někdy se pod fundamentálními diagramy rozumí i závislosti q(k), q(v) či jiné. Obecně lze říct, že rychlost toku musí s hustotou klesat. Nejjednodušším možným modelem je lineární závislost mezi rychlostí a hustotou. Takový model formuloval v roce 1935 Greenshields, který navrhl závislost   k , v = vf 1 − kj kde vf je rychlost na volné komunikaci a kj hustota dopravy při kongesci. Podle tohoto vztahu se vozidla na volné komunikaci pohybují rychlostí vf . Se zvětšující se hustotou k rychlost vozidel klesá. Jakmile hustota dosáhne hodnoty kj , rychlost vozidel klesá na nulu a formuje se kongesce. Hustota toku je v tomto modelu dána vztahem   k q = k v = k vf 1 − . kj Fundamentální diagramy pro Greenshieldsův model jsou znázorněny na obrázku 1.1. U obrázku 1.1 je maximální tok dosažen pro hodnotu k = kj /2[1, str. 18,34].

10

2. Mikroskopické modely Mikroskopické modely se, na rozdíl od makroskopických, zabývají vlastnostmi souvisejícími s každým jednotlivým vozidlem. Zatímco makroskopické modely jsou inspirovány mechanikou kontinua, mikroskopické modely mají svůj protějšek v molekulové fyzice, ve které se chování systému jako celku vysvětluje jako výsledek složité interakce jednotlivých konstituentů systému.V dopravě to znamená, že výsledný dopravní tok je výsledkem chování každého jednoho řidiče. Interakce mezi molekulami mohou být značně složité, ale fyzikálně lze předpokládat, že každá daná molekula neinteraguje se všemi ostatními molekulami, ale pouze s těmi, které se nacházejí v jejím bezprostředním okolí. Analogicky je chování řidiče ovlivněno nikoli všemi vozidly nacházejícími se na komunikaci, nýbrž pouze vozidly, které má na dohled. Cílem mikroskopických modelů je co nejvěrněji popsat chování řidiče pomocí obyčejných diferenciálních rovnic. Parametry, na které se řidič zaměřuje, jsou rychlost, zrychlení, vzdálenost mezi vozidly a funkce těchto proměnných (např. čas do kolize). Řidič si tyto parametry sjednotí a vyhodnotí podle svých vlastních zkušeností i dalších okolností (například parametrů svého vozidla). Při ovládání vozidla pak může řidič provádět jednotlivé úkony s lehkostí, obratností, a to v závislosti na svých předchozích zkušenostech. Již v 50. letech 20. století se stal proces porovnávání řidičových „vstupůÿ a „výstupůÿ aktuální. Problémem ovšem bylo, že nebyla nalezena žádná adekvátní přenosová funkce, která by těmto vstupům a výstupům odpovídala. Díky tomu vznikl nápad vytvořit modely, které nastavují zrychlení jednotlivých vozidel v závislosti na jejich okolí, tzv. Car Following Models (CFM). Při využití mikroskopických modelů se dá lépe simulovat chování řidiče ve stresu, předjíždění a dalších situací, při kterých se změní parametry pouze jediného vozidla[8, str. 160,161]. První modely CFM se zabývaly rozestupem a časovou vzdáleností mezi vozidly. Jejich název proto zní Safety distance models 1 . V těchto modelech se vychází z předpokladu, že řidič dodržuje bezpečnou vzdálenost k vozidlu před sebou. Pro tyto modely je základem klasický Newtonův pohybový zákon. V některých modelech se dokonce počítá tato bezpečná vzdálenost pouze jako vzdálenost nutná k zabránění nehodě pokud přední vozidlo náhle sníží svoji rychlost. Pipesův model vycházel z předpokladu, že bezpečná mezera mezi vozidly závisí lineárně na rychlosti, konkrétně, že se zvýší o délku vozidla při nárůstu rychlosti druhého vozidla oproti prvnímu o 16 km/h. V roce 1981 byl testován další model, Gippsův. V něm má řidič vozidla zajištěno, že se nesrazí s vozidlem před ním, pokud je časová vzdálenost mezi oběma vozidly větší než 1,5 násobek reakčního času řidiče a pokud řidič správně odhadl rychlost brždění vozidla před sebou[8, str. 164-166]. Jiný náhled na věc zaujala skupina výzkumníků ze společnosti General Motors (GM), která se zabývala CFM modely na principu akce a reakce. Existuje pět modelů GM, podle toho, jak se postupně zdokonalovala měření a matematické výpočty. Základním předpokladem všech těchto modelů je rovnice odezvy2 Odezva = f (citlivost, podnět).

(2.1)

Odezva je v GM modelech vždy nahrazena zrychlením aktuálního vozidla zatímco podnět je reprezentován rozdílem rychlostí mezi aktuálním a následujícím vozidlem. Rozdíly mezi jednotlivými GM modely spočívají v určení citlivostní funkce. Tato myšlenka pochází z klasické newtonovské mechaniky, kde je zrychlení částice považováno za její reakci k podnětu způsobenou silou, která je obdržena interakcí s jinými částicemi v systému. 1 Model

založený na bezpečných vzdálenostech mezi vozidly equation

2 Stimul-response

11

V prvním GM modelu byl citlivostní člen α považován za konstantní. Pohybová rovnice tu má tvar v˙ n (t + ∆t) = α [vn−1 (t) − vn (t)] ,

(2.2)

kde ∆t je reakční čas. Záleží tedy pouze na rozdílu rychlostí vedoucího a aktuálního auta. Podnětový člen může být kladný, když vedoucí auto jede rychleji vn−1 > vn , což vyústí v zrychlování aktuálního vozidla. Podnětový člen bude nulový za podmínky vn−1 = vn , takže zrychlení aktuálního vozidla bude nulové a jeho rychlost konstantní. Za podmínky vn−1 < vn bude aktuální vozidlo zpomalovat. Experimentálně bylo dokázáno, že citlivostní člen je konstantní pouze pro úzké spektrum reakčních časů a bezpečných mezer, proto byl výzkumníky model vylepšen z rovnice (2.2) na tvar v˙ n (t + ∆t) = α1 [vn−1 (t) − vn (t)]

nebo

v˙ n (t + ∆t) = α2 [vn−1 (t) − vn (t)] ,

(2.3)

v závislosti na bezpečnostní mezeře. Pokud byla vozidla u sebe velmi blízko, použila se vysoká hodnota citlivostního členu, tzn. α1 , pokud od sebe byla vozidla daleko, použila se nízká hodnota α2 . Záhy se zjistilo, že je velký problém určit, kdy už jsou od sebe vozidla daleko a kdy jsou ještě blízko. Byly proto prováděny další pokusy, jejichž výsledkem byl experimentální vztah pro citlivostní člen α α0 α0 = , α0 ∈ h0, 1i , (2.4) α= d xn−1 (t) − xn (t) kde d je vzdálenost mezi vozidly, α0 je parametr závislý na d a α. Čím blíže tedy vozidla k sobě budou, tím vyšší hodnotu bude mít citlivostní člen. Pokud se tento vztah dosadí do (2.3), vychází v˙ n (t + ∆t) =

α0 [vn−1 (t) − vn (t)] , xn−1 (t) − xn (t)

(2.5)

což je zároveň rovnice třetího GM modelu. Citlivostní člen je závislý na bezpečnostní vzdálenosti a na parametru α0 . Jeho rozměr je m/s. Ve čtvrtém modelu GM se opět vylepšuje citlivostní člen přidáním rychlosti následujícího vozidla. Tento model vychází z toho, že pokud se zrychluje dopravní proud, řidič bude mnohem citlivější na rozdíl rychlostí. Rovnice toho modelu vypadá takto v˙ n (t + ∆t) =

α0 [vn (t + ∆t)] [vn−1 (t) − vn (t)] , xn−1 (t) − xn (t)

(2.6)

citlivostní člen se skládá ze tří komponent, rychlosti vozidla, parametru α0 (který je bezrozměrný) a bezpečnostní mezery. Poslední, pátý, model se snaží vylepšit a zobecnit citlivostní člen, tím, že jak čitatel, tak jmenovatel jsou umocněny na určité koeficienty, které se mění v závislosti na typu komunikace. Rovnice pátého modelu má tvar v˙ n (t + ∆t) =

α0 [vn (t + ∆t)]m l

[xn−1 (t) − xn (t)]

[vn−1 (t) − vn (t)] ,

(2.7)

[8, str. 167-173]. Modely GM předpokládají, že řidič reaguje na relativně malé změny při své jízdě, dále také reaguje na vozidlo před ním i tehdy, pokud už je od něj velmi vzdálené. a řidič také nebude 12

∆x

z´ ona bez reakce

z´ ona reakce

z´ ona reakce

0

∆v

Obrázek 2.1: Zóny pro reakci řidičů v závislosti na rozdílu rychlostí a poloh reagovat při stejných rychlostech vozidel. Tyto nedostatky mohou být odstraněny buďto rozšířením GM modelů novými režimy - jízda na volné silnici, pohotovostní brzdění anebo použitím psycho-fyziologického modelu, což je nový přístup k CFM, který byl představen v 70. letech. V psycho-fyziologických modelech jsou použity určité prahy, při jejichž dosažení změní řidič druhého vozidla své chování. Prahy jsou v nejjednodušším případě definovány dva, a to vzájemná rychlost a vzdálenost vozidel. Z obrázku 2.1 vyplývá, že když rozdíl rychlostí obou vozidel dosáhne určitého prahu, předpokládá určitou reakci ze strany druhého vozidla. Tato reakce ovšem nenastane, pokud je vzájemná vzdálenost vozidel taková, že také dosáhne určité hranice. Tyto modely jsou implementovány v známých mikrosimulačních softwarech jako je VisSim, který používá Wiedermannův model či AIMSUN [1, str. 25]. Posledním typem CFM modelů jsou modely s fuzzy logikou, která kvantifikuje výroky, například výrok „Jede moc rychle.ÿ Pokud tento výrok splní koeficient příslušnosti, tak se aktivuje rychlé brzdění. V předchozích modelech řidiči museli znát aktuální rychlost a polohu ostatních řidičů či další parametry, což zde odpadá. V těchto modelech se pouze kontroluje, zda automobil splnil určitý výrok a pokud ano, tak se vykonají jednotlivé příkazy. Vozidlo tedy může být v několika množinách zároveň a tím pádem se vykonává více příkazů. To je hlavní výhodou tohoto přístupu. Byly snahy fuzifikovat GM modely, to se ovšem do této doby nepodařilo[7].

13

3. Dynamické systémy Z matematického hlediska je studium mikroskopických dopravních modelů založeno na teorii obyčejných diferenciálních rovnic. Vhodný teoretický rámec pro tuto analýzu je teorie dynamických systémů, což jsou soustavy obyčejných diferenciálních rovnic prvního řádu. V této kapitole uvedeme některé základní poznatky této teorie. Systém se nazývá stabilní tehdy, když se při malé změně počátečních podmínek se chování systému v budoucích stavech změní pouze málo. Naopak nestabilita znamená, že sebemenší změna počátečních podmínek způsobí velký rozdíl v jeho dalším chování. Ukazuje se, že velmi důležitou součástí matematických systémů nejsou pouze jejich rovnice, ale též počáteční podmínky. Význam této vlastnosti systémů rozpoznal James Clerk Maxwell1 , který uvedl, že nepatrná nejistota ve znalosti současného stavu znemožní přesnou předpověď stavu budoucího. Z toho plyne, že aby se daly předpovídat budoucí stavy systému, musely by být známy počáteční podmínky s dokonalou přesností. V tu chvíli by byl popis deterministického systému k užitku. To se ovšem stát nemůže[2, str. 78,79]. S tímto názorem se později ztotožňoval i francouzský matematický fyzik Henri Poincaré, který studoval problém, jak klasický deterministický popis může vést k pravděpodobnostnímu vyjádření. Avšak vznik kvantové mechaniky2 , která byla schopna odpovědět na některé otázky týkající se chaosu, přibrzdil další studium citlivosti na počátečních podmínkách. Navíc v té době ještě nebyl dostatečný matematický aparát tyto úlohy řešit. To přišlo až s příchodem výpočetní techniky, kdy bylo možné řešit rozsáhlé numerické simulace. Dalším matematikem, který se zajímal o determinismus byl Pierre Simon de Laplace. Ten stanovil deterministický model světa, který umožňuje na základě počátečních podmínek a přírodních zákonů stanovit všechny podstatné jevy budoucnosti i minulosti. Poincaré ovšem jeho tvrzení vyvrátil tím, že minimální změny přírodních zákonů mohou vést k úplně nepředvídatelným výsledkům. Některé nelineární dynamické systémy vykazují deterministický chaos. To znamená, že ačkoliv je systém popsán rovnicemi velice přesně a neobsahuje nahodilé parametry, v důsledku citlivosti na počáteční podmínky pak působí chaoticky. Hovoří se o tzv. „Motýlím efektuÿ. Jak později ukázal Lorenz, chaotičnost se objevuje už i v poměrně jednoduchých nelineárních systémech. Zjistil, že v jeho zjednodušeném modelu atmosféry existují chaotické trajektorie, které dnes nazýváme Lorenzovým atraktorem, viz obrázek 3.1 [6, str. 7].

3.1

Planární dynamický systém

V této části bude pro jednoduchost uvažován autonomní dynamický planární systém, který je popsán pouze dvěmi proměnnými x1 a x2 . Poté se rovnice dynamického systému (4.1) zjednoduší na x˙ 1 = f1 (x1 , x2 )

x˙ 2 = f2 (x1 , x2 ).

(3.1)

Fázovým prostorem pro ně pak bude rovina. Pro vyšetřování stability je klíčové nalézt kritické body. Říká se jim také pevné, stacionární, rovnovážné či singulární. Pokud budou počáteční 1 „Je-li stav věcí takový, že nekonečně malá variace přítomného stavu změní pouze o nekonečně malou veličinu stav v nějakém budoucím čase, pak podmínky systému, ať už je v klidu či v pohybu, nazveme stabilními; povede-li však nekonečně malá změna přítomného stavu ke konečnému rozdílu ve stavu systému v konečném čase, nazveme podmínky systému nestabilní.[2, str. 158]ÿ 2 kvantová mechanika je doposud jediná teorie, která ukazuje, že svět je náhodný na úrovni mikrosvěta; příroda ve skutečnosti není deterministická

14

20

0

-20

40

30

20

10

0 -10 0 10 20

Obrázek 3.1: Lorenzův atraktor. podmínky systému nastaveny do kritického bodu, nic se nestane, systém zůstane ve stejné poloze. Podle toho, jak bude systém reagovat na nepatrnou změnu počátečních podmínek, se rozlišují typy kritických bodů, viz dále. K nalezení kritických (pevných) bodů se položí jednotlivé funkce rovny nule f1 (x1c , x2c ) = f2 (x1c , x2c ) = 0, kde x1c a x2c jsou souřadnice kritického bodu[10, str. 156]. Ke zjištění stability jednotlivých kritických bodů se zkoumá chování systému v jejich okolí. Za tímto účelem zavedeme nové proměnné  a δ vztahy x1 = x1c + δ, x2 = x2c + ,

δ  1,

(3.2)

  1.

Toto vyjádření se dosadí do rovnice (3.1) a zanedbají se všechny členy vyšších řádů v δ a , čímž se dostavá linearizovaný systém. Zde platí HartmanGrobmanův teorém, který říká, že chování dynamického systému v blízkosti kritického bodu bude kvalitativně stejné i poté, co se provede jeho linearizace. V okolí kritických bodů se použije zjednodušený zápis systému ve tvaru δ˙ = a δ + b ,

(3.3)

˙ = c δ + d ,

15

Λ1 < 0, Λ2 < 0

4

4

2

2

0

0

y

y

Λ1 > 0, Λ2 > 0

-2

-2

-4

-4

-4

-2

0

2

4

-4

x

-2

0

2

4

x

Obrázek 3.2: Trajektorie pro λ1,2 > 0

Obrázek 3.3: Trajektorie pro λ1,2 < 0

kde koeficienty a,b,c a d jsou parciálními derivacemi v Jacobiho matici plynoucí z Taylorova rozvoje   ∂f1 ∂f1    ∂x ∂x2  1 a b   J= = (3.4) . c d  ∂f ∂f2  2 ∂x1 ∂x2

3.2

Klasifikace kritických bodů

Jednotlivé kritické body lze klasifikovat na základě vlastních čísel λ1 a λ2 matice (3.4), která ovlivňují systém v okolí kritického bodu. Možnosti jsou následující: • reálné kořeny λ1 a λ2 – λ1,2 > 0 – nestabilní uzel, trajektorie vycházejí z počátku a směřují do nekonečna (obrázek 3.2) – λ1,2 < 0 – stabilní uzel, trajektorie směřují do počátku, (obrázek 3.3) – λ1 · λ2 < 0 – sedlový bod, trajektorie jsou přitahovány k ose s kladnou hodnotou a odpouzovány od osy se zápornou hodnotou λ (obrázek 3.4) • komplexné sdružené kořeny λ1 a λ2 (λ1,2 = α ± iβ) – α = 0 – střed, trajektorie jsou kružnice se středem v počátku (obrázek 3.5 a 3.6) – α < 0 – stabilní ohnisko, trajektorie jsou přitahovány k počátku (obrázek 3.7) – α > 0 – nestabilní ohnisko, trajektorie směřují do nekonečna (obrázek 3.8)

16

Λ1 < 0, Λ2 > 0

4

y

2

0

-2

-4

-4

0

-2

2

4

x

Obrázek 3.4: Trajektorie pro λ1 · λ2 < 0

Α = 0, Β < 0

Α = 0, Β > 0

4

4

2

2

0

0

-2

-2

-4

-4

-4

-2

0

2

4

-4

Obrázek 3.5: Trajektorie pro α = 0, β < 0

-2

0

2

4

Obrázek 3.6: Trajektorie pro α = 0, β > 0

17

Α < 0, Β > 0

Α > 0, Β > 0

4

4

2

2

0

0

-2

-2

-4

-4

-4

-2

0

2

4

-4

Obrázek 3.7: Trajektorie pro α < 0, β > 0

-2

0

2

4

Obrázek 3.8: Trajektorie pro α > 0, β > 0

Pokud se vlastní čísla rovnají, λ1 = λ2 , pak se kritický bod nazývá degenerovaný, v opačném případě se nazývá singulární. Vlastní čísla λ1,2 splňují rovnici J · e = λ · e, kde každému vlastnímu číslu λi náleží vlastní vektor ei . Geometrický význam vlastního vektoru spočívá v tom, že představuje invariantní systém, tedy směr, který se působením matice J nemění. Invariantní směr (vlastní vektor matice J) se nazývá varieta. Ta je stabilní pro λ < 0, kdy se trajektorie přibližuje k počátku, nestabilní varieta odpovídá kladné hodnotě λ[10, str. 160-173].

18

4. Numerické metody Často se uvádí, že numerické metody jsou současně vědou i uměním. Využívají se hlavně tehdy, kdy není možné najít řešení určitého problému analyticky nebo řešení lze nalézt, ale je to buď časově nebo ekonomicky náročné nebo není nutné, aby řešení bylo úplně přesné. V našem světě jsme si zvolili elementární funkce (například polynomy, goniometrické, logaritmus), které ovšem zdaleka nedokáží popsat všechna možná řešení různých problémů. Proto v první řadě využíváme numerické metody na přibližné určení správného řešení. Jiným problémem stále zůstavá, jak přijít ke správným diferenciálním rovnicím pro popis toho určitého systému.

4.1

Dynamický systém

Uvažujme systém (obecně nelineárních) diferenciálních rovnic prvního řádu ve tvaru x˙ 1 = f1 (x1 , x2 , . . . , xm , t) x˙ 2 = f2 (x1 , x2 , . . . , xm , t) .. .

(4.1)

x˙ m = fm (x1 , x2 , . . . , xm , t), pro které chceme nalézt řešení numericky[14, str. 124]. Tyto rovnice jsou zapsány v tzv. úplném tvaru. Rovnice také můžeme zapsat v kratším vektorovém tvaru x˙ = f (x, t),

(4.2)

kde x = (x1 , x2 , . . . , xm ). Podobně pak f = (f1 , f2 , . . . , fm ). Složky vektoru x resp. f označíme xa ,

fa ,

kde a = 1 . . . m.

(4.3) 1

Soustavu (4.1) nebo (4.2) budeme nazývat dynamickým systémem . Takovýto samotný systém není pro výpočet postačující, protože pokud nejsou známy doplňující neboli počáteční podmínky. V tu chvíli není jasné odkud má daná numerická metoda vycházet.2 Pro výpočet proto bude také nutno znát m doplňujících podmínek x1 (0) = x01 ,

x2 (0) = x02 ,

...

, xm (0) = x0m .

(4.4)

Horní index u proměnné bude označovat její hodnotu v příslušném čase. V tomto případě značí hodnotu v čase t = 0. Může se také stát, že je zadaná diferenciální rovnice m-tého řádu y (m) = f (y (m−1) , y (m−2) , . . . , y, t),

(4.5)

Pro řešení takové rovnice se zavedou nové nezávislé proměnné: y1 = y, y2 = y˙ 1 , .. .

(4.6)

ym = y (m−1) , 1 či

Cauchyho systémem vozidlo, které nemá zadanou svou počáteční rychlost a dráhu.

2 Například

19

a pomocí těchto proměnných se problém převede na systém diferenciálních rovnic (4.1), který již lze numericky řešit y˙ 1 = y2 , y˙ 2 = y3 , .. .

(4.7)

y˙ m = f (y1 , y2 , . . . , ym , t), [14, str. 124-127]. Tento převod se používá například při zjednodušení diferenciální rovnice 2. řádu a = s¨, kde a je zrychlení a s je dráha. Přidáním nové nezávislé proměnné v, tj. rychlosti, vzniknou 2 diferenciální rovnice 1. řádu a = v, ˙

v = s. ˙

Toto zjednodušení neboli redukce řádu je pouze formální a z praktického hlediska neulehčuje řešení diferenciální rovnice. Pod „formálnímÿ zjednodušením se myslí to, že získaná soustava rovnic prvního řádu je ve skutečnosti stejně komplikovaná jako výchozí diferenciální rovnice. Důvod, proč se preferuje tato redukce, spočívá v tom, že dynamický systém má názornou geometrickou interpretaci a všechny použité numerické metody jsou sestaveny pro rovnice prvního řádu.

4.2

Eulerova metoda

Nejjednodušší numerickou metodou je metoda Eulerova. Je sice nejméně přesná, ale také nejjednodušší a lze na ní vysvětlit základní principy využívané i u složitějších metod. První věcí je diskretizace času, tj. zavedení posloupnosti časových okamžiků tk = hk,

(4.8)

kde h značí tzv. integrační krok či diferenci. U Eulerovy metody je tato diference konstantní. Princip Eulerovy metody vychází z definice derivace x˙ a = lim

h→0

xa (t + h) − xa (t) = fa (x(t), t). h

Nyní se pouze nahradí spojitý čas jeho diskretizovanou verzí xk+1 − xka a = fa (x(tk ), tk ) h a po úpravě vychází jeden iterační krok Eulerovy metody xk+1 = xa (tk ) + h fa (x(tk ), tk ) ≡ xka + h fa (xk , tk ). a

(4.9)

Důležité je si uvědomit, že tento iterační krok je absolutně nezávislý na konkrétní funkci fa . Eulerův integrátor dostane na začátku informace o počátečních podmínkách x0a a o funkci fa a pouze na základě těchto informací dokáže nalézt řešení v libovolném čase. Iterační schéma

20

Euler krok

x

x

x2 x3 x1

~τ1

~τ0

y = f (x)

x2 x3 x1

~τ0

~τ1 ~τ

~τ1

x = f (t) RK2 krok

t1

t2

t3

t

t1 t2 t3 t b) metoda Runge-Kutta 2. řádu

a) Eulerova metoda

Obrázek 4.1: Grafické znázornění numerických metod zůstavá tedy pro jakýkoliv problém stejné, jediné co ve schématu (4.9) na studovaném problému závisí, je fa (x(tk ), tk ). Na první pohled je jasné, že čím víc se bude zmenšovat integrační krok h, tím přesnější se dostane výsledek. U numerických metod se rozlišují dva základní druhy chyb. Jak je vidět na obrázku 4.1a, Eulerova metoda pracuje tak, že v bodě [t1 , x1 ] nalezneme tečnu a pohneme s ní ve směru vektoru τ 0 do bodu [t2 , x2 ] a v bodě [t2 , x2 ] opět nalezneme tečnu a pohneme s ní ve směru vektoru τ 0 . Rozdíl t2 − t1 resp. t3 − t2 je právě diference h(4.8). Tímto způsobem se ale postupně vzdaluje od správného řešení. Eulerova metoda je tedy přesná pouze tehdy, když její integrační krok se limitně blíží nule, h → 0, [9, str. 132,133].

4.3

Metoda Runge-Kutta

Runge-Kutta pracuje poněkud přesněji oproti Eulerově metodě. U Eulerovy metody vypočteme souřadnice bodu [t2 , x2 ] (viz. obrázek 4.1a) tak, že t2 = t1 + h

x2 = x1 + y(t ˙ 1 ) = x1 + f (x0 )h.

Metoda Runge-Kutta se snaží tento problém „obejítÿ. Stejně jako v Eulerově metodě se nejprve posuneme z bodu [t1 , x1 ] ve směru tečny ~τ0 , ale pouze s polovičním integračním krokem h t2 = t1 +

h 2

h x2 = x1 + f (x0 ) . 2

V tomto bodě se opět vypočítá derivace funkce, což je směrnice ~τ1 . A tato tečna se nyní zpátky posune do bodu [t1 , x1 ]. Nyní se udělá „průměr tečenÿ. ~τ =

~τ1 + ~τ2 2

V tomto směru se potom Runge-Kutta posune o celý integrační krok h. Jak je vidět z obrázku 4.1b, řešení metodou Runge-Kutta druhého řádu je přesnější než základní Eulerovou metodou. Přesnějších výsledků lze tedy dosáhnout nejen snižováním diference h, ale také vypočtením hodnot derivací v bodech, kde se nepřičítá celá hodnota diference, ale jen její zlomek. U metody Runge Kutta druhého řádu se zapíše základní iterační krok takto [9, str. 145] 21

xk+1 = xka + a

1 (k1 + k2 ) , 2

kde  k2 = hfa

xka

k1 h + , tk + 2 2



a k1 = h fa (xka , tk ). Tento proces je graficky znázorněn na obrázku 4.1b. Metody Runge-Kutta se používají nejen 2. řádu, ale i víceřádové. Jejich výhodou bezpochyby je zvýšení přesnosti a to tím, že se vypočítávají derivace ve vícero bodech, které jsou rozprostřeny na intervalu htk , hi. Nejčastější ze všech těchto metod je metoda 4. řádu. Její iterační krok vypadá následovně[9, str. 146]: 1 xk+1 = xka + h (k1 + 2k2 + 2k3 + k4 ) , a 6

(4.10)

k1 = fa (xk , tk ),  1 k2 = fa xk + k1 , tk + 2  1 k3 = fa xk + k2 , tk + 2
k k4 = fa x + k3 , tk .

(4.11)

kde

 h , 2  h , 2

Existují i další metody 4. řádu [9, str. 145-148], zde ovšem plně postačí schéma 4.10. U metod vyšších řádů výpočetní náročnost algoritmu převažuje nad zvýšením přesnosti.

22

5. Intelligent Driver Model Jeden z představitelů „car following modelsÿ je také Intelligent Driver Model (IDM), model inteligentního řidiče. Tento model patří do skupiny mikroskopických modelů, ve kterých se sleduje každý účastník na komunikaci zvlášť a ne komunikaci jako celek, nicméně nepatří do žádné ze skupin uvedených v sekci 2. Dále patří do tzv. deterministických modelů, to znamená, že při opakujících se stejných vstupech bude systém vydávat stejné výstupy.1 Tento model uvažuje pouze jednoproudou komunikaci (bez možnosti předjíždění)[13]. Mezi základní vstupy do IDM se zařazuje u každého účastníka (auta) jeho rychlost a poloha. Mezi parametry tohoto modelu se uvádí bezpečný rozestup mezi vozidly, průměrné zrychlení vozidla, brzdová konstanta, preferovaná rychlost vozidla a bezpečná doba zastavení před jiným vozidlem. Některé tyto parametry můžou mít jednotlivá auta jiné (např. preferovaná rychlost). Jako výstupy se potom dostávají závislosti polohy, rychlosti či zrychlení jednotlivých aut na čase[5, str.2]. Rozhodování, zda účastník bude v danou chvíli zrychlovat či zpomalovat, závisí pouze na jeho rychlosti a pozici a také na rychlosti a pozici účastníka před ním. V každém čase se potom dají sledovat jednotlivé vlastnosti každého účastníka (rychlost, poloha. .). Mezi hlavní výhody IDM modelu patří[5, str.2,3]: • je do velké míry nekolizní • všechny parametry modelu mají „rozumnouÿ interpretaci, jsou empiricky měřitelné a jsou známy rozsahy jejich hodnot • z fundamentálního diagramu a vlastností stability modelu lze snadno a samostatně kalibrovat data do empirických dat • hodí se pro numerické simulace • akcelerace se může lišit od decelerace, např. v případě nehody auto začne brzdit razantněji (asymetrie modelu) • patří mezi tzv. akcelerační modely, které mají výhodu v tom, že hodnota aktuálního zrychlení nebude nesmyslná hodnota při dobře nastavených parametrech simulace • je vhodný pro popis tří hlavních typů proudů na silnicích – volný proud – účastník provozu jede úplně sám a není ničím omezován, tohoto proudu lze dosáhnout např. na závodním okruhu – plynulý proud – účastník provozu jede plynule, tj. musí se sice přizpůsobovat ostatním účastníkům, ale tito účastníci mají podobné či stejné parametry jako on sám (nízká intenzita i hustota dopravy) – dopravní zácpa – účastník provozu je bržděn účastníkem, který jede před ním (nízká intenzita a vysoká hustota dopravy) 1 Opakem deterministických modelů jsou modely stochastické, kdy do pohybových rovnic vstupují i náhodné veličiny simulující nepředvídatelné prvky chování řidičů či jiné náhodné vlivy; u těchto modelů jsou výsledky nereprodukovatelné

23

5.1

Matematický popis

Pro výpočet zrychlení v IDM je použit vztah v˙ n = afree + aint ,

(5.1)

kde se aktuální zrychlení účastníka vypočítá jako součet zrychlení při volném proudu afree a zrychlení, ve kterém je zahrnut i účastník jedoucí před ním aint . Pro volný proud je tedy zrychlení aint = 0. Zrychlení při volném proudu je definováno vztahem " δ #  vn , (5.2) afree = a 1 − vnfree kde a je maximální zrychlení, δ je exponent zrychlení, vn je aktuální rychlost vozidla a vnfree je rychlost, které chce vozidlo dosáhnout. Pokud v rovnici (5.2) bude zlomek vn /vfree roven jedné, auto se již dostalo na svou preferovanou rychlost vfree a celkové zrychlení afree bude rovno nule. Naopak pokud bude vn rovno nule, znamená to, že afree bude přímo rovno a, a auto pojede svým maximálním zrychlením. Pokud vn bude menší než vfree , pak platí že, čím vyšší hodnotu bude mít nastavenou parametr δ, tím se dosáhne vyšší hodnota afree . Pokud vn bude větší než vfree , bude platit pravý opak. Auto jedoucí při volném proudu je ale nereálná představa, proto se do rovnice (5.1) přidává interakční člen, který přímo závisí na parametrech vozidla jedoucí před naším účastníkem. Tento člen je roven  ∗ 2 s , (5.3) aint = −a n sn kde a je maximální zrychlení, sn = xn−1 − xn je rozdíl poloh vozidel a s∗n s∗n = s0 + vn T +

vn (vn − vn−1 ) √ 2 ab

(5.4)

je preferovaný odstup. s0 je součet minimální bezpečné vzdálenosti mezi vozidly a délkou následujícího vozidla, T je bezpečnostní časový odstup mezi vozidly a b je konstanta brzdění. Situace je schématicky znázorněna na obrázku 5.1. Člen úměrný vn − vn−1 je kvadratický v proměnné vn a obsahuje součin vn vn−1 . Tyto nelineární členy popisují interakci, v tomto případě dvou vozidel. Pokud tedy bude vzdálenost sn mezi auty vysoká, potom zlomek s∗n /sn bude mít velmi malou hodnotu a interakční člen aint lze zanedbat, takže vozidlo se pohybuje jako ve volném proudu. Pokud rychlost auta vn bude vyšší než rychlost auta vn−1 , pak interakční člen s∗n bude kladný a požadovaný odstup se zvětší, takže zrychlení auta v˙ n se musí snížit. Hodnoty, které se pro simulaci běžně používají jsou uvedeny v tabulce 5.1.

5.2

Model MOBIL

Při simulaci vozidel je také potřeba brát v úvahu, že vedoucí vozidlo může jet nižší rychlostí než vozidla za ním nebo se na vozovce může nacházet nějaká překážka. Z tohoto důvodu byl také vymyšlen předjížděcí model MOBIL[13]. Pro řidiče je předjíždění na komunikacích i za normálních okolností považováno za problém. Řidič se vždy rozhoduje na základě několika faktorů, nejdůležitějším z nich je bezpečnost. Dalším 24

vn n−1

n sn

vn−1

s∗ Obrázek 5.1: Proměnné jednotlivých aut Parametr Rychlost, které chce auto dosáhnout vfree Časový odstup T Minimální vzdálenost s0 Průměrné zrychlení a Konstanta brzdění b

Hodnoty pro automobil 120 km/h 1.5 s 2.0 m 0.3 m/s2 3.0 m/s2

Hodnoty pro kamion 80 km/h 1.7 s 2.0 m 0.3 m/s2 2.0 m/s2

Tabulka 5.1: Hodnoty používané pro simulaci Zdroj:[13] z nich jsou například jeho vlastní zkušenosti. Dále se pak řidič rozhoduje podle výhodnosti vedlejšího pruhu, což by mělo vést k primárně zvýšení jeho rychlosti a tím i rychlejšímu dosažení jeho preferované rychlosti. Některá vozidla by měla preferovat pravý jízdní pruh i tehdy, když je pro ně výhodnější jet v levém pruhu2 , jiná vozidla by mohla volit zase levý jízdní pruh. Všechny tyto faktory se snaží model MOBIL zachytit.

5.2.1

Popis modelu

Podstatné pro předjíždění jsou dvě kritéria, která musí být splněna: • druhý jízdní pruh musí být pro automobil více atraktivní, tj. „motivační kritériumÿ • změna jízdního pruhu je bezpečná, nehrozí nebezpečí kolize, tj. „bezpečnostní kritériumÿ Pro lepší názornost je situace znázorněna na obrázku 5.2. Jednotlivé obdelníky značí automobily, první index značí pozici automobilu, druhý index pak jízdní pruh. Vozidlo, které bude uvažovat o změně pruhu, je označeno písmenem B. Bezpečnostní kritérium Bezpečnostní kritérium se stará o to, aby nedošlo k situaci, že vozidlo přejede do vedlejšího pruhu a vozidlo, které v tomto pruhu již jede, bude muset prudce brzdit, případně ani zabrzdit nedokáže. Závisí tudíž na automobilu, které jede ve vedlejším pruhu za ním (E). U vozidla E se tedy spočte teoretické zrychlení an+1,2 podle vztahu (5.1) při situaci, kdyby se auto B nacházelo v levém pruhu. Pokud by hodnota zrychlení byla menší než hodnota předem stanovená konstanta, vozidlo E by nebylo schopno zareagovat na přejetí vozidla B a mohlo by dojít i ke srážce. Toto zrychlení tedy musí být větší než bezpečnostní konstanta bsafe . Vyjádřeno matematicky an+1,2 > −bsafe .

(5.5)

2 Jako příklad se dá uvést kamion, který jede rychleji než kamion před ním, ovšem řidič kamionu ví, že bude následovat stoupání, předjíždel by pomalu a proto zpomalí a zůstane ve stávajícím pruhu.

25

an+1,2

an−1,2

E

D

an+1,1

an,1

A

an−1,1

B

C

Obrázek 5.2: Schéma situace při výpočtu předjíždění Pokud tedy tato nerovnice platí, je splněno bezpečnostní kritérium. Motivační kritérium Bezpečná změna pruhu pro vozidlo ještě nutně nemusí znamenat změnu pruhu. Od toho musí být splněno motivační kritérium. Toto kritérium vychází z předpokladu, že vozidlo chce dosáhnout co nejrychleji své preferované rychlosti. Z počátku se může zdát, že toto kritérium je zde proto, aby umožnilo plynulou jízdu vozidla, tj. aby se zrychlení vozidla měnilo co nejméně. To ovšem není totéž. Pro motivační kritérium je zapotřebí, aby byla splněna nerovnice (5.6). Na levé straně se vcelku logicky nachází rozdíl zrychlení automobilu ve vedlejším pruhu a zrychlení automobilu ve stávajícím pruhu. Pro větší přesnost se většinou na levé straně ještě uvádí konstanta br , která určuje preferenci pruhu. Na dvouproudové silnici by vozidla měla jet vždy napravo a měnit jízdní pruh pouze pokud chtějí předjíždět, a ovšem hned poté by se měla zařadit zpátky. V modelu MOBIL se toho docílí nastavením konstanty br na příslušnou hodnotu, která bude pro pravý jízdní pruh kladná a pro levý záporná. Na pravé straně nerovnice by čistě teoreticky stačila hodnota athr , která eliminuje vliv „bláznivé jízdyÿ, tj. jízdy, kdy by auto neustále přejíždělo z jednoho pruhu do druhého. Nachází se zde ale také rozdíl zrychlení pro vozidlo E. Počítá se nejprve zrychlení vozidla E v závislosti na automobilu před ním a také na automobilu B, pokud by změnilo pruh. To, jak moc se řidič vozidla B ohlíží na vozidlo E, je určeno „zdvořilostníÿ proměnnou p. Pokud tedy bude hodnota p rovna nule, je jasné, že vozidlo B se vůbec neohlíží na vozidlo E a přejede do druhého pruhu. Takovýto řidič se dá nazvat pirátem silnic. Jediný náhled vozidla B na vozidlo E by tedy zůstal v bezpečnostním kritériu. Naopak pokud je hodnota p rovna jedné znamená to, že řidič se v podstatě neustále dívá do zpětného zrcátka a bude předjíždět pouze tehdy, pokud tím v podstatě vůbec nezpomalí vozidlo E. Pokud je hodnota p vyšší, vozidlo nemění svůj jízdní pruh či mění jen pokud je to nezbytně nutné. Normální hodnota na silnici pro p je 0,5. Ale pokud je například na silnici nějaká překážka tak se hodnota p snižuje někdy až k nule. Na dálnici se hodnota pohybuje od 0 do 0,5. Pokud se položí parametry p = 1 a athr = 0, změna jízdního pruhu se vykoná přesně v případě, že by se předjetím pouze zvýšilo zrychlení, a žádné auto by nezpomalilo. Takto také vznikl název pro tento model MOBIL Minimizing Overall Braking decelerations Induced by Lane changes3 . 3 Minimalizace

celkového zpomalení vyvolané změně jízdního pruhu.

26

Matematicky se podmínka motivačního kritéria zapíše jako br + an,2 − an,1 > p ∗ [an+1,2 (D) − an+1,2 (B)] + athr .

27

(5.6)

Část II

Vypracování

28

6. Popis programu Program má v sobě implementovaný model IDM a MOBIL, dále je zde možnost nastavení cyklické dálnice. Jako výstup programu se vytvoří soubor, který obsahuje čas, polohu a rychlost jednotlivých vozidel, případně dalších vlastností, které se dají donastavit. Základní okno programu je vidět na obrázku 6.1. To je rozděleno na několik částí, z nichž každá má jinou funkci. JControlPanel informuje u vybraného vozidla o jeho vlastnostech a o vlastnostech vozidel před a za vozidlem v tom samém pruhu. CarSelector slouží k vybrání libovolného auta. To se stane aktivním jak v JControlPanelu, tak i v panelu JHighway, kde se zobrazí doprostřed. JHighway pak zobrazuje samotnou dálnici, kde je v této situaci vidět aktuální vybrané vozidlo 0 (značené žlutou barvou), ostatní vozidla jsou označena bílou barvou. Horní posuvník v JHighway mění měřítko dálnice, na začátku je nastaveno na 40m, max. hodnota může dosáhnout až 500m.

6.1

Převod kartézských souřadnic

Prvním problémem, který je v programu řešen, je převod souřadnic. Pro jednoznačné určení každého bodu sítě dálnice se použije lineární kombinace dvou vektorů u = (ux , uy ) a v = (vx , vy ). Pro nějaký bod P [x, y] tedy bude platit, že jeho souřadnice se dají napsat ve tvaru x = α ux + β v x ,

(6.1)

y = α uy + β vy ,

kde α a β jsou koeficienty, určující přesnou pozici bodu v souřadném systému. V tuto chvíli je nutné si uvědomit, že počítač nemá stejný souřadný systém jako je kartézský. Počítač bere jako počátek souřadnic bod Op a od něj dopočítává všechny ostatní body. Proto souřadnice bodu O budou po přepočtu O[borderX, h − borderY ]. Oblast dálnice je určena kosodelníkem jehož vrcholy jsou ABDO. Jednotlivé složky vektorů se vypočítají jako rozdíly souřadnic bodů A a O pro vektor ~u a bodů B a O pro vektor ~v . Tím se docílí toho, že odteď se budou jednotlivé body počítat od bodu O a ne od bodu Op . Tyto vektory je nutné vydělit vzdáleností, která se bude lišit v závislosti na pozici a měřítku. U osy y se měřítko měnit nebude a bude to šířka dálnice. U osy x je problém poněkud složitější. Zde se mění měřítko i pozice, od které se bude daný úsek dálnice sledovat. Ax − Ox xa − x0 Bx − Ox vx = yb − y0

ux =

Ay − Oy xa − x0 By − Oy vy = yb − y0

uy =

(6.2)

Pokud by se souřadný systém měnil tak, že by osa x a osa y na sebe byly vzájemně kolmé, hodnota uy a vx by se potom rovnala nule. To znamená, že by se rovnice (6.1) zjednodušily tak, že nová hodnota souřadnice by závisela pouze na jednom vektoru. Tím pádem by se daly transformovat souřadnice x a y zvlášť. V případě izometrické dálnice na sebe osy kolmé nejsou a nulová je pouze hodnota uy . Je tedy vždy nutné transformovat obě souřadnice zároveň. To, co se v případě změny měřítka bude měnit je hodnota xa a x0 . Hodnota x0 nastavuje pozici, od které se bude dálnice sledovat a hodnota xa se bude rovnat konci sledované oblasti. Rozdíl tedy určuje vzdálenost, která bude zobrazovaná. U hodnot yb a y0 zústavá hodnota stejná. Rozdíl určuje šířku dálnice. 29

Obrázek 6.1: Snímek obrazovky programu Při počítání vykreslení bodu P se nejdříve musí zkontrolovat, zda souřadnice [x,y] leží v zobrazované oblasti. Pokud ano, dopočtou se hodnoty α a β, které jsou hodnoty rozdílu x − x0 resp. y − y0 . Nové souřadnice pak budou mít hodnotu vypočtenou pomocí (6.1). Všechny body a souřadnice jsou znázorněny na obrázku 6.2.

6.2

Balík Utilities

V tomto balíku jsou dvě třídy. První z nich Conversions obsahuje funkce primárně pro převod jednotek, druhá Vectors pak funkce pro práci s vektory.

6.2.1

Třída Conversions

Třída Conversions je pomocná třída pro převod jednotek rychlosti, hodnoty úhlů. Hodnotu, ve kterém pruhu se nachází automobil, je uložena ve formě čísla, proto třída obsahuje funkci pro převod tohoto čísla na text. Seznam funkcí: • public static double Radians(int degrees) - převod stupňů na radiány • public static double KmHtoMS(double kmh) - převod rychlosti z km/h na m/s • public static double MStoKmH(double ms) - převod rychlosti z m/s na km/h • public static String LaneToStr(int lane) - v závislosti na pruhu funkce vrátí textový řetězec

30

Op

w D

B

P [x,y] h ~v = (vx , vy )

O[borderX, h − borderY ]

borderX

[x0 , y0 ]

~u = (ux , uy )

A [xa , y0 ]

borderY

Obrázek 6.2: Popis bodu pomocí lineární kombinace

6.2.2

Třída Vectors

V programu budou automobily jezdit po izometrické dálnici. Pro tento účel je vytvořena třída Vectors, která zajišťuje základní práci s vektory. Java totiž tyto funkce primárně neobsahuje. Seznam funkcí: • public static double[] Vector(double[] A, double[] B) - ze dvou bodů (vstupních parametrů) se vypočte vektor • public static double[] Multiply(double[] v, double k) - každá složka vektoru v bude vynásobena konstantou k • public static double[] LinComb( double alpha, double[] u, double beta, double[] v) - vypočítání vektoru, který je lineární kombinací vektoru u a v ~a = α ~u + β ~v • public static double[] LinComb( double[] r0, double alpha, double[] u, double beta, double[] v) - stejná funkce jako předešlá, s tím rozdílem, že se dá nastavit vzdálenost od počátku ~a = r0 + α ~u + β ~v

31

6.3

Balík mathcomps

6.3.1

Třída JDrawPanel

Základem pro třídu JDrawPanel bude klasický JPanel z třídy Swing1 . Všechny funkce tohoto panelu se týkají vykreslení izometrické dálnice. Dále jsou zde funkce ke kreslení čar, tak i vyplnění určité oblasti barvou. Všechny souřadnice budou zadávány v kartézských souřadnicích a v této třídě budou převedeny na počítačové (viz kapitola 6.1). K nakreslení přerušované čáry slouží funkce fillRhomboid. Seznam funkcí: • public JDrawPanel () - konstruktor, nastavuje barvu pozadí panelu na černou • public void setScale (double[] Origin, double[]A, double[]B, double x0, double y0, double xA, double yB) - nastavení měřítka panelu. Tato funkce pracuje na základě výpočtu vektorů u a v ve směrech určenými body A a B. Poté se tyto vektory vynásobí koeficientem závisejícím na stanoveném měřítku. Každý bod kosodelníku je pak jednoznačně určen lineární kombinací vektorů u a v. • public int[] transform(double x, double y) - převod souřadnic bodů z kartézských do počítačových. Použije se k tomu funkce z balíku Vectors LinComb, kde se nové (počítačové) souřadnice bodu spočtou jako lineární kombinace vektoru u a v. Koeficienty α a β jsou vzdálenosti od počátku souřadnic (bod[x0 , y0 ]). • public void move(Graphics g, double x, double y) - uloží do proměnné lastPoint pozici odkud se posléze bude kreslit čára • public void line(Graphics g, double x, double y) - nakreslí čáru z bodu lastPoint do bodu [x,y] • public void fillRhomboid( Graphics g, double a1, double b1, double a2, double b2, double a3, double b3, double a4, double b4) - vyplní oblast určenou vstupními parametry.

6.4

Třída ControlPanel

Třída ControlPanel je opět rozšířením třídy JPanel. Tentokrát je ovšem vytvořen panel, který poskytuje informace o jednotlivých vozidlech. • public void createGUI() - funkce pro vykreslení všech popisků, nastavení pozadí • public void arrangeGUI() - nastavení rozmístění prvků na ControlPanelu • public void refreshData() - funkce, která se volá při každé změně, vypisuje u vybraného auta rychlost, pozici a preferovanou rychlost, dále pak informace rychlost a vzdálenost k autu následujícímu a předchozímu • public void setCar(int i) - volá se při výběru vozidla, poté se ihned zavolá refreshData() • public void initialize(carsList) - inicializace celého panelu 1 Tato

třída je určena pro vytváření grafických komponent a prvků ve formuláři.

32

6.5

Třída CarList

Tato třída je základním kamenem programu. Obsahuje v sobě třídu CarInfo, která nastavuje parametry každému jednotlivému automobilu, které je vygenerováno. Všechny parametry jsou stejné jako v kapitole 5.1. Jelikož jsou všechny hodnoty proměnných (tz. vlastnosti automobilu) již nastavené, pro vytvoření nového auta stačí použít příkaz 1

Carinfo car = new CarInfo(defaultImageFile);

a v proměnné car již jsou uložena všechna data o automobilu. Např. v proměnné 1

car.lane

je uložena hodnota o aktuálním pruhu vozidla. Třída CarInfo • double x, double v, double vFree, double a, double b, double L, double T, int lane - proměnné určující parametry vozidla, jejich vysvětlení viz. 5.1 • double aThreshold, double p, double bSafe - parametry předjízděcího modelu MOBIL, 5.2 • boolean overtaking - proměnná určující zda je povoleno předjíždění • public CarInfo(String imageFile) - vytvoří novou instanci vozidla s obrázkem, který je parametrem konstruktoru • public int nextLane() - vrací druhý jízdní pruh než ve kterém jede vozidlo Třída CarList • public CarInfo[] cars - pole aut, které se při zavolání konstruktoru CarList naplní vozidly, a pracuje se s ním v celé třídě • double hwLength - určuje délku dálnice, její hodnota se využije pouze tehdy, když je nastavena cyklická dálnice • boolean cyclic - určuje zda je nastavena cyklická dálnice, tj. pokud vozidlo dojede na konec dálnice, která je určena hodnotou hwLength, tak pokud je cyclic nastavena na true, vozidlo se objeví na začátku dálnice • public CarList(int numberOfCars) - konstruktor, který vytvoří pole automobilů dle vstupní hodnoty numberOfCars • public int getCount () - vrací délku pole cars • public double sqr (double x) - vrací druhou mocninu čísla x • public double power(double x, int n) - umocní číslo x na n • public int findLeader(int i, int lane) - funkce vrací hodnotu indexu vozidla, které jede před vybraným (indexem i ) v závislosti na jízdním pruhu (vstupní parametr lane) • public int findFollower(int i, int lane) - totožná funkce jako předchozí, s tím rozdílem, že vrací index vozidla jedoucí za vybraným Tyto funkce pracují na principu vypočítávání polohy auta s indexem i a porovnávání s polohou všech ostatních vozidel. Pokud je hodnota menší než u předchozího vozidla v poli, 33

uloží se rozdíl poloh do pomocné proměnné a stejně tak jako index tohoto vozidla, jinak se nestane nic. Rozdíl mezi FindLeader a FindFollower je pouze v tom, jak se vypočítává vzdálenost mezi těmito vozidly. V případě funkce FindLeader je to poloha ostatních automobilů a od ní odečtena poloha automobilu s indexem i. U funkce FindFollwer je tomu právě naopak. Posledním problémem je to pokud automobil pojede první, resp. poslední. Index vozidla je proto na začátku funkce nastaven na hodnotu -1, a pokud se po prozkoumání poloh ostatních aut zjistí, že žádná vzdálenost není kladná, index zůstane nastaven na -1, což značí, že automobil jede první resp. poslední. • public double getFrontDistance(int i, int j) - funkce počítající vzdálenost nejbližšího vozidla jedoucího před aktuálním • public double getBackDistance(int i, int j) - funkce počítající vzdálenost nejbližšího vozidla jedoucího za aktuálním Tyto funkce jsou nutné hlavně pokud je nastavena cyklická dálnice. Jinak by při zobrazování vzdáleností stačilo odečíst souřadnice obou vozidel. Pokud ovšem vozidlo jede těsně před koncem dálnice a druhé hned na začátku, tak by rozdíl vzdáleností nefungoval, resp. vycházel by nesmysl. • public double acceleration (int i, int j) - funkce vypočítávající přesné zrychlení aint podle modelu IDM a rovnic uvedených v podkapitole 5.1. Parametry jsou indexy aut, mezi kterými probíhá interakce. (Pokud auto jede jako první, zrychlení bude nulové, aint = 0.) • public void nextStep() - hlavní funkce, vykonávající všechny potřebné úkony během jednoho posunu Eulerovy metody. Pro každé vozidlo se najde vozidlo před ním, za ním, ve vedlejším pruhu, spočte se jeho zrychlení, dále se ověřuje to, zda automobil bude měnit svůj jízdní pruh. Poté se změní jeho rychlost a dráha, případně i jízdní pruh.

6.6

Třída Highway

Tato třída se stará o grafické uživatelské rozhraní celého programu. Pokud se třída Highway rozšíří třídou JFrame, která bude mít implementované události na určitou změnu, získá se prvotní vzhled okna. Třída JHighWay, která je rozšířením JDrawPanelu, se stará o správné vykreslování dálnice. Pokud se okno spouští poprvé nebo se mění rozměry okna, zavolá se vytvořená funkce arrangeGUI, která přepočítá umístění jednotlivých komponent ve formuláři. Třída JHighWay • double hwLength, double markerLength, double markerWidth, double hwWidth, int borderX, int borderY - konstanty určující vykreslení dálnice, jejich význam je znázorněn na obrázku 6.3 • double position - důležitá proměnná určující pozici, od které se začne dálnice vykreslovat • double zoom - proměnná, jejíž hodnota stanovuje počet zobrazovaných metrů (měřítko dálnice) • private void drawHighway(Graphics g) - tato funkce se stará o vykreslení dálnice včetně přerušované čáry uprostřed • public void paint(Graphics g) - funkce volá funkci DrawHighway a poté už se stará o vykreslení aut a popisků, které se právě nacházejí v zobrazované oblasti dálnice

34

Třída HighWay • public Highway() - konstruktor, který zobrazuje formulář, provádí základní nastavení okna (minimalizace, správné ukončení) a volá funkci pro vykreslení obsahu formuláře • private void createGUI() - funkce pro vytvoření objektů ve formuláři. Formulář je rozložen na několik částí. První z nich je panel, na který se vykresluje dálnice, nad ním je posuvník určující měřítko dálnice. Další část tvoří panel, ve kterém se zobrazují informace o právě vybraném automobilu a informace o předchozím a následujícím automobilu, poslední jsou pak zbylá tlačítka, která zajišťují generování, spouštění a přerušování simulace. • private void arrangeGUI() - nastavuje rozměry jednotlivým komponentám ve formuláři. • public void componentResized(ComponentEvent e); public void componentShown (ComponentEvent e) - při zobrazení či při změně rozměrů okna se automaticky zavolá funkce arrangeGUI • public void run() - pokud běží vlákno, provede se nextStep, ControlPanel obnoví data, a pomocí funkcí třídy FileCreator se zapíšou informace o automobilech do souboru • public double equilibriumVelocity(double vLow,double vHigh,int i, double distance) - funkce, která pracuje úplně stejně jako metoda sečen. V tomto případě hledá rychlost v rozsahu vLow až vHigh, při které bude zrychlení nulové. • public double y(double x,int i, double distance) - funkce která pro dané x nalezne funkční hodnotu y Reakce na komponenty umístěné na formuláři • GenerateButton - po stisknutí tohoto tlačítka se vygenerují jednotlivá vozidla, nastaví se jim aktuální a preferovaná rychlost a startovní pozice, inicializuje se ControlPanel z informacemi a vytvoří se soubor pro zápis (viz. 6.7) • CarSelected - nastaví se auto vybrané ze seznamu jako aktuální a podle něj se překleslí dálnice • ZoomChange - při změně měřítka se pouze mění hodnota proměnné zoom • DrawingButton - při stisku se přestane vykreslovat dálnice, simulace ovšem běží dál • RunButton - při startu simulace se vytvoří nové vlákno a spustí se; pokud simulace běží, vlákno se ukončí

6.7

Třída FileCreator

Tato třída slouží pro výpis informací do souboru. • public void initialize(CarList carList, String gnuFileName) - vytvoření a nastavení souboru pro externí program • public void finalize() - funkce starající se o řádné ukončení práce se souborem • public void writeStatus(CarList carList) - zápis času a pozice každého automobilu z pole carList

35

w [x3 , y3 ]

[x2 , y2 ]

h markerLength

[x0 , y0 ]

[x1 , y1 ]

borderX

borderY

Obrázek 6.3: Konstanty použité pro vykreslení dálnice

36

7. Zdrojový kód V této příloze uvádíme kompletní výpis zdrojového kódu programu popsaného v této bakalářské práci.

7.1

Hlavní program Listing 7.1: Highway.java

1 2

import mathcomps.JDrawPanel; import utilities.Conversions;

3 4 5 6 7

import import import import

javax.swing.*; java.awt.*; java.awt.event.*; java.util.Random;

8 9

public class Highway extends JFrame implements ComponentListener, Runnable {

10 11 12

private class JHighWay extends JDrawPanel {

13 14 15 16 17

double hwLength = 1000; /* in meters */ double position = 0; /* position of the left edge */ double zoom = 40; /* meters per screen */ double viewAngle = Conversions.Radians(45); /* in degrees */

18 19 20

double markerLength = 1.5; /* marker between lanes */ double markerWidth = 1;

21 22 23 24

double hwWidth = 10; /* width in meters */ int borderX = 50; /* border along the road in pixels */ int borderY = 70;

25 26

int w, h;

27 28 29

Polygon road; Font numbersFont = new Font("Times New Roman", Font.PLAIN, 20);

30 31 32 33 34 35 36 37 38 39 40 41 42

private void drawHighway(Graphics g) { g.setColor( Color.GRAY ); /* corners of the road */ int[] xs = new int[4]; int[] ys = new int[4]; xs[0] = borderX; ys[0] = h - borderY; //left down xs[1] = w - 2*borderX; ys[1] = h - borderY; //right down xs[2] = w - borderX; ys[2] = borderY; //right up xs[3] = 2*borderX; ys[3] = borderY; road = new Polygon(xs, ys, 4); g.fillPolygon(road); g.setClip(road);

37

double[] O = {xs[0], ys[0]}; double[] A = {xs[1], ys[1]}; double[] B = {xs[3], ys[3]};

43 44 45 46

setScale(O, A, B, position, -hwWidth/2, position+zoom, hwWidth/2); g.setColor(Color.WHITE);

47 48 49

/* markers */

50 51

int im = (int)(position/(2*markerLength) - 1); double xm; do { xm = (2*im+1) * markerLength; fillRhomboid(g, xm, -markerWidth / 2, xm + markerLength, -markerWidth / 2, xm + markerLength, markerWidth / 2, xm, markerWidth / 2 ); im++; } while (xm < position + zoom);

52 53 54 55 56 57 58 59 60 61 62 63 64

}

65 66 67 68

public void update(Graphics g) { paint(g); }

69 70 71

public void paint (Graphics g) { super.paint(g);

72 73

w = getWidth(); h = getHeight();

74 75 76 77 78 79 80

position = (carList==null)?0:carList.cars[selectedCar].x-zoom/2; drawHighway(g); move(g,position, -hwWidth/2); line(g,position, hwWidth/2); if (carList==null) return;

81 82 83 84 85 86 87 88 89

for (int i=0; i < carList.getCount(); i++ ) { CarList.CarInfo car = carList.cars[i]; if ((car.x>=position)&(car.x<=position+zoom)) { int[] coords = transform(car.x, 0); g.drawImage(car.image, coords[0], (car.lane==1)?coords[1]-30:coords[1]-90, null); } }

90 91 92 93 94 95

g.setFont(numbersFont); g.setClip(null); Rectangle roadBounds = road.getBounds(); for (int i=0; i < carList.getCount(); i++) if ((carList.cars[i].x >= position)&(carList.cars[i].x<=position+zoom)) {

38

int[] coords = transform(carList.cars[i].x, 0); g.setColor((i==selectedCar)?Color.YELLOW:Color.WHITE); g.drawString(""+i, coords[0] + 30, roadBounds.y + roadBounds.height + 30);

96 97 98

}

99 100 101

} }

102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121

public void run() { while (animator!=null) { carList.nextStep(); if (drawing) highWay.repaint(); controlPanel.refreshData(); timeLabel.setText("Time: t = " + Math.round(carList.time)); fileWriteSteps++; if (fileWriteSteps>fileWriteStepsMax) { fileWriteSteps = 0; fileCreator.writeStatus(carList); } try { // Thread.sleep(80); } catch (Exception e) { System.out.println("Thread problem"); } } } Thread animator = null;

122 123 124 125 126 127 128 129 130 131 132

JHighWay highWay; JPanel contentPane; JButton runBtn, generateBtn, drawingBtn, specialBtn; CarList carList = null; JComboBox carSelector; JScrollBar zoomBar; JLabel zoomLabel; JLabel timeLabel; ControlPanel controlPanel; FileCreator fileCreator;

133 134 135 136

int fileWriteStepsMax = 100; int fileWriteSteps = fileWriteStepsMax; boolean drawing = true;

137 138 139

int selectedCar;

140 141 142 143 144 145 146 147 148

public Highway() { super("Highway simulation"); setExtendedState( getExtendedState()|JFrame.MAXIMIZED_BOTH ); setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE); setBounds(100, 100, 400, 400 ); createGUI(); setVisible(true); // loading picture

39

149 150

}

151 152 153 154

public static void main (String[] args) { new Highway(); }

155 156 157 158 159 160 161

private void createGUI() { contentPane = new JPanel(); contentPane.setBackground( Color.DARK_GRAY ); contentPane.setLayout( null ); addComponentListener( this );

162 163 164 165

setContentPane( contentPane ); highWay = new JHighWay(); add(highWay);

166 167 168 169 170

runBtn = new JButton("Start simulation"); runBtn.addActionListener(new RunButton()); runBtn.setEnabled(false); add(runBtn);

171 172 173 174

generateBtn = new JButton("Generate cars"); generateBtn.addActionListener(new GenerateButton()); add(generateBtn);

175 176 177 178 179 180

carSelector = new JComboBox(); carSelector.setEditable(false); carSelector.setEnabled(false); carSelector.addActionListener(new CarSelected()); add(carSelector);

181 182 183 184 185 186 187

zoomBar = new JScrollBar(JScrollBar.HORIZONTAL, 40, 5, 20, 500); zoomBar.addAdjustmentListener(new ZoomChange()); add(zoomBar); zoomLabel = new JLabel("Meters / screen = " + highWay.zoom); zoomLabel.setForeground(Color.YELLOW); add(zoomLabel);

188 189 190 191

controlPanel = new ControlPanel(); controlPanel.setBackground(Color.DARK_GRAY); add(controlPanel);

192 193 194 195

timeLabel = new JLabel("Time: t = "); timeLabel.setForeground(Color.YELLOW); add(timeLabel);

196 197 198 199

drawingBtn = new JButton("Stop drawing"); drawingBtn.addActionListener(new DrawingButton()); add(drawingBtn);

200 201

specialBtn = new JButton("");

40

specialBtn.addActionListener(new SpecialButton()); add(specialBtn);

202 203 204

}

205 206

private void arrangeGUI() { zoomBar.setBounds(50, 30, getWidth() - 100, 20); zoomLabel.setBounds(50, 10, getWidth() - 100, 20); highWay.setBounds(50, 50, getWidth()-100, getHeight()/3); runBtn.setBounds(50, getHeight()/3 + 70, 200, 30); generateBtn.setBounds(300, getHeight()/3 + 70, 200, 30); drawingBtn.setBounds(520, getHeight()/3 + 70, 200, 30); specialBtn.setBounds(750, getHeight()/3 + 70, 200, 30); carSelector.setBounds(250, getHeight()/3 + 130, 100, 30); controlPanel.setBounds(50, getHeight()/3 + 170, getWidth()/3, (int)(1.5*getHeight()/3)); controlPanel.arrangeGUI(); timeLabel.setBounds(50 + (getWidth()-100)/2, 10, getWidth()-100, 20); // specialBtn.setText(SpecialEventsList[specials]); }

207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223

public void componentHidden(ComponentEvent e) {} public void componentResized(ComponentEvent e) { arrangeGUI(); } public void componentMoved(ComponentEvent e) {} public void componentShown(ComponentEvent e) { arrangeGUI(); }

224 225 226 227 228 229 230 231 232

private class RunButton implements ActionListener { public void actionPerformed (ActionEvent e){ if (animator==null) { runBtn.setText("Stop simulation"); fileCreator.initialize(carList, "cars"); animator = new Thread(Highway.this); animator.start(); } else { animator.interrupt(); animator = null; runBtn.setText("Start simulation"); fileCreator.finalize();

233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248

}

249

}

250 251

}

252 253 254

private class GenerateButton implements ActionListener { public void actionPerformed (ActionEvent e) {

41

int N = 10; carList = new CarList(N);

255 256 257

double double double double

258 259 260 261

distance = 100; vFree = Conversions.KmHtoMS(120); vLast=0; xLast=0;

262

//equidistantly distributed for (int j=0;j < N-1; j++){ CarList.CarInfo car = carList.cars[j]; car.x = j*distance; car.vFree = vFree; vLast = car.v = equilibriumVelocity(0, 2*vFree, 0, distance); }

263 264 265 266 267 268 269 270

// we set up the leading car separately CarList.CarInfo car = carList.cars[N-1]; car.x = (N-1)*distance; car.v = car.vFree = vLast;

271 272 273 274 275 276

fileCreator = new FileCreator();

277 278

carSelector.removeAllItems(); for (int i=0; i < carList.getCount(); i++) carSelector.addItem("Car " + i); carSelector.setEnabled(true );

279 280 281 282 283

highWay.repaint(); runBtn.setEnabled(true); controlPanel.initialize(carList); controlPanel.arrangeGUI(); controlPanel.setCar(0); generateBtn.setEnabled(false);

284 285 286 287 288 289

}

290 291

}

292 293 294 295 296 297 298 299

private class CarSelected implements ActionListener { public void actionPerformed (ActionEvent e ) { selectedCar = carSelector.getSelectedIndex(); controlPanel.setCar(selectedCar); highWay.repaint(); } }

300 301 302 303 304 305 306 307

private class ZoomChange implements AdjustmentListener { public void adjustmentValueChanged (AdjustmentEvent e) { highWay.zoom = e.getValue(); zoomLabel.setText("Meters / screen = " + e.getValue()); if ((animator!=null)&(!animator.isAlive())) { highWay.repaint(); }

42

}

308 309

}

310 311 312 313 314 315 316 317 318 319 320 321 322

private class DrawingButton implements ActionListener { public void actionPerformed(ActionEvent e) { if (drawing) { drawingBtn.setText("Restore drawing"); drawing = false; } else { drawingBtn.setText("Stop drawing"); drawing = true; } } }

323 324 325 326 327

int specials = 0; public static final String[] SpecialEventsList = { "Slow down leading car", "Accelerate the same car"};

328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350

private class SpecialButton implements ActionListener { public void actionPerformed(ActionEvent e) { CarList.CarInfo car; switch (specials) { case 0: /* leading car slows down */ car = carList.cars[carList.getCount()-1]; car.vFree = Conversions.KmHtoMS(50); break; case 1: /* accelerate it back */ car = carList.cars[carList.getCount()-1]; car.vFree = Conversions.KmHtoMS(250); car.a = 1; car.p = 0; car.rightBias=0; car.leftUnBias=0; break; } specials = (specials==1)?0:specials+1; specialBtn.setText(SpecialEventsList[specials]); } }

351

public double y(double x,int i, double distance){ CarList.CarInfo car = carList.cars[i]; double afree = car.a * (1 - carList.power(x/car.vFree, car.delta));

352 353 354 355

double s = distance; double sStar = car.L + x*car.T;// return afree - car.a * carList.sqr(sStar / s);

356 357 358 359 360

}

43

361

public double equilibriumVelocity(double vLow,double vHigh,int i, double distance){ double nextIteration; double n=0; double epsilon = 1E-10; do{ n++; nextIteration = vHigh((y(vHigh,i,distance)*(vHigh-vLow))/ (y(vHigh,i,distance)-y(vLow,i,distance))); if ((y(nextIteration,i,distance)*y(vLow,i,distance))<0){ vHigh = nextIteration; } else { vLow = nextIteration; } } while(Math.abs(y(nextIteration,i,distance))>=epsilon); System.out.println("Method of secants converged after " + n + " steps"); return nextIteration; }

362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383

}

Listing 7.2: ControlPanel.java 1

import utilities.Conversions;

2 3 4

import javax.swing.*; import java.awt.*;

5 6

public class ControlPanel extends JPanel {

7 8 9 10 11 12 13 14 15 16 17

private private private private private private private private private private

JLabel JLabel JLabel JLabel JLabel JLabel JLabel JLabel JLabel JLabel

currentCarLabel; speedLabel; positionLabel; vFreeLabel; nextCarLabel; nextCarSpeedLabel; nextCarDistanceLabel; prevCarLabel; prevCarSpeedLabel; prevCarDistanceLabel;

private private private private

int carIndex = 0; CarList carList; boolean initialized = false; int width, height;

18 19 20 21 22 23 24 25 26 27 28

public void initialize(CarList carsList) { createGUI(); // carList = carsList;

44

initialized = true;

29 30

}

31 32 33 34 35 36

public void setCar(int i) { carIndex = i; // refreshData(); }

37 38 39 40 41 42 43 44

public void createGUI() { setBackground(Color.BLACK); setLayout(null); // currentCarLabel = new JLabel("Selected car:"); currentCarLabel.setForeground(Color.YELLOW); add(currentCarLabel);

45

positionLabel = new JLabel("Position: x = "); add(positionLabel); speedLabel = new JLabel("Speed: v = "); add(speedLabel); vFreeLabel = new JLabel("Free speed: vf = "); add(vFreeLabel);

46 47 48 49 50 51 52

nextCarLabel = new JLabel("Next car: "); nextCarLabel.setForeground(Color.YELLOW); add(nextCarLabel);

53 54 55 56

nextCarDistanceLabel = new JLabel("Distance: s = "); add(nextCarDistanceLabel);

57 58 59

nextCarSpeedLabel = new JLabel("Speed: v = "); add(nextCarSpeedLabel);

60 61 62

prevCarLabel = new JLabel("Previous car: "); prevCarLabel.setForeground(Color.YELLOW); add(prevCarLabel);

63 64 65 66

prevCarDistanceLabel = new JLabel("Distance: s = "); add(prevCarDistanceLabel);

67 68 69

prevCarSpeedLabel = new JLabel("Speed: v = "); add(prevCarSpeedLabel);

70 71 72 73

}

74 75 76 77 78 79 80 81

public void arrangeGUI() { if (!initialized) return; width = getWidth(); height = getHeight(); currentCarLabel.setBounds(20, 0, width-40, 30); positionLabel.setBounds(20, 30, width-40, 30);

45

speedLabel.setBounds(20, 60, width-40, 30); vFreeLabel.setBounds(20, 90, width-40, 30); nextCarLabel.setBounds(20, 120, width-40, 30); nextCarDistanceLabel.setBounds(20, 150, width-40, 30); nextCarSpeedLabel.setBounds(20, 180, width-40, 30); prevCarLabel.setBounds(20, 210, width-40, 30); prevCarDistanceLabel.setBounds(20, 240, width-40, 30); prevCarSpeedLabel.setBounds(20, 270, width-40, 30);

82 83 84 85 86 87 88 89 90

}

91 92 93 94 95 96 97 98 99 100 101 102

public void refreshData() { if (!initialized) return; CarList.CarInfo car = carList.cars[carIndex]; currentCarLabel.setText("Selected car: " + carIndex); positionLabel.setText("Position x = " + Math.round(car.x) + " m"); speedLabel.setText("Speed v = " + Math.round(Conversions.MStoKmH(car.v)) + " km/h"); vFreeLabel.setText("Free speed vf = " + Math.round(Conversions.MStoKmH(car.vFree)) + " km/h"); int iNext = carList.findLeader(carIndex, car.lane);

103 104 105 106 107 108 109 110 111 112 113 114 115 116

if (iNext < 0) { nextCarLabel.setText("this car is the leading car in the " + Conversions.LaneToStr(car.lane)); nextCarDistanceLabel.setText(""); nextCarSpeedLabel.setText(""); } else { CarList.CarInfo nextCar = carList.cars[iNext]; nextCarLabel.setText("Next car: " + iNext); nextCarDistanceLabel.setText("Distance: s = " + Math.round(nextCar.x - car.x) + " m"); nextCarSpeedLabel.setText("Speed: v = " + Math.round(Conversions.MStoKmH(nextCar.v)) + " km/h"); }

117 118

int iPrev = carList.findFollower(carIndex, car.lane);

119 120 121 122 123 124 125 126 127 128 129 130 131 132

if (iPrev < 0) { prevCarLabel.setText("this car is the last car in the " + Conversions.LaneToStr(car.lane)); prevCarDistanceLabel.setText(""); prevCarSpeedLabel.setText(""); } else { CarList.CarInfo prevCar = carList.cars[iPrev]; prevCarLabel.setText("Previous car: " + iPrev); prevCarDistanceLabel.setText("Distance: s = " + Math.round(prevCar.x - car.x) + " m"); prevCarSpeedLabel.setText("Speed: v = " + Math.round(Conversions.MStoKmH(prevCar.v)) + " km/h"); }

133 134

if ((iNext>=0)&(carList.cars[iNext].x - car.x<0))

46

nextCarDistanceLabel.setText("Distance: s = " + Math.round(carList.hwLength+carList.cars[iNext].x - car.x) + " m"); if ((iPrev>=0)&(carList.cars[iPrev].x-car.x>0)) prevCarDistanceLabel.setText("Distance: s = " + Math.round(-(carList.hwLength-carList.cars[iPrev].x + car.x)) + " m");

135 136 137 138 139 140

}

141 142

}

Listing 7.3: CarList.java 1 2 3

import javax.imageio.ImageIO; import java.awt.*; import java.io.File;

4 5

public class CarList {

6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

public class CarInfo { double x = 0; double v = 0; double vFree = 0; double a = 0.3; double b = 3; double L = 7; double T = 1.5; double p = 0.5; double rightBias = 0.3; double leftUnBias = -0.1; int lane = 1; double aThreshold = 0.2; double bSafe = 4; int delta = 4; public Image image; public String imageName;

24 25 26 27 28 29 30 31 32 33

public CarInfo(String imageFile) { try { imageName = imageFile; image = ImageIO.read(new File(imageFile)); } catch (Exception e) { System.out.println("Error reading " + imageFile); } }

34 35 36 37

public int nextLane() { return (lane==1)?2:1; } }

38 39 40 41 42 43

private int n; double time = 0; private String defaultImageFile = "car.png"; public CarInfo[] cars; double hwLength = 1000;

47

44 45 46 47 48 49 50 51

boolean cyclic = true; public CarList(int numberOfCars) { n = numberOfCars; cars = new CarInfo[n]; for (int i = 0; i
52 53 54 55

public int getCount () { return cars.length; }

56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71

private CarInfo copyCar(CarInfo source) { CarInfo target = new CarInfo(source.imageName); target.x = source.x; target.v = source.v; target.vFree = source.vFree; target.a = source.a; target.b = source.b; target.L = source.L; target.T = source.T; target.lane = source.lane; target.aThreshold = source.aThreshold; target.bSafe = source.bSafe; target.p = source.p; return target; }

72 73 74 75 76 77 78 79

private CarList copyCars(CarList source) { int n = source.getCount(); CarList target = new CarList(n); for (int i =0; i
80 81

double dt = 0.01;

82 83 84 85

public double sqr (double x) { return x*x; }

86 87 88 89 90 91 92

public double power(double x, int n) { double y = 1; for (int i = 0; i < n; i++) y *= x; return y; }

93 94 95 96

public int findLeader(int i, int lane) { if (i<0) return -1;

48

97 98 99 100 101 102 103 104 105 106 107 108 109 110

double minDist = -1; //negative distance corresponds to leading car double dist; int lead = -1; for (int j = 0; j < getCount(); j++) if ((i!=j)&(cars[j].lane==lane)) { dist = cars[j].x - cars[i].x; if ((dist<0)&cyclic) dist = (hwLength-cars[i].x) + (cars[j].x - 0); if ((dist>0)&((minDist<0)|(dist<minDist))) { minDist = dist; lead = j; } } return lead; }

111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130

public int findFollower(int i, int lane) { if (i<0) return -1; double minDist = -1; //negative distance corresponds to the last car double dist; int follow = -1; for (int j = 0; j < getCount(); j++) if ((i!=j)&(cars[j].lane==lane)) { dist = cars[i].x - cars[j].x; if ((dist<0)&(cyclic)) dist = (hwLength-cars[j].x) + (cars[i].x-0); if ((dist>0)&((minDist<0)|(dist<minDist))) { minDist = dist; follow = j; } } return follow; }

131 132 133 134 135 136 137 138 139 140 141 142 143

public double acceleration (int i, int j) { if ((j < 0)|(i<0)) return 0; CarInfo car = cars[i]; double s = cars[j].x car.x>0?cars[j].x - car.x:hwLength-car.x+cars[j].x; double sStar = car.L + car.v*car.T + car.v * (car.v - cars[j].v) / (2 * Math.sqrt(car.a * car.b)); return - car.a * sqr(sStar / s); }

144 145

boolean overtaking = true;

146 147 148 149

public void nextStep() { CarList old = copyCars(this); for (int i = 0; i < old.getCount(); i++) {

49

CarInfo car = cars[i]; int leader = findLeader(i, car.lane); //leading car in the current lane int leaderNext = findLeader(i, car.nextLane()); //leading car in the next lane int follower = findFollower(i, car.nextLane()); //following car in the next lane double aFree = (car.vFree==0)?0:car.a * (1 - power(car.v/car.vFree, car.delta)); double acc1 = old.acceleration(i, leader); //acceleration in the current lane double acc2 = old.acceleration(i, leaderNext); //acceleration in the next lane // acceleration of previous car in the next lane double accB = old.acceleration(follower, findLeader(follower, car.nextLane())); // hypothetical acceleration of previous car if the current car // changes the lane double accB2 = old.acceleration(follower, i); // we prefer to ride in the right lane double right = (car.lane==2)?car.rightBias:car.leftUnBias;

150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171

boolean incentive = right + acc2 - acc1 > car.p*(accB - accB2) + car.aThreshold; boolean politeness = accB2 > - car.bSafe;

172 173 174 175

double acc;

176 177

if (overtaking&incentive&politeness) { acc = acc2; car.lane = car.nextLane(); } else { acc = acc1; } car.v += (aFree + acc) * dt; car.x += car.v * dt; if (cyclic) if (car.x>=hwLength){ car.x = car.x-hwLength; } time += dt;

178 179 180 181 182 183 184 185 186 187 188 189 190

}

191

}

192 193

}

Listing 7.4: FileCreator.java 1 2

import java.io.*; import java.util.Scanner;

3 4

public class FileCreator {

5 6 7

private String gnuFile; private PrintWriter out;

50

8

private PrintWriter densityFile;

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

public void initialize(CarList carList, String gnuFileName) { try { gnuFile = gnuFileName; FileWriter outFilePositions = new FileWriter(gnuFileName + ".pos"); PrintWriter outPositions = new PrintWriter(outFilePositions); outPositions.println("reset"); outPositions.println("set style data lines"); outPositions.print("plot "); for (int i = 0; i < carList.getCount(); i++) { outPositions.print("’" + gnuFileName + "’ using 1:" + (i+2) + " notitle"); if (i < carList.getCount() - 1) { outPositions.print(", "); } } outPositions.close(); // FileWriter outFile = new FileWriter(gnuFileName); out = new PrintWriter(outFile); } catch (Exception e) { System.out.println("Error creating " + gnuFile); }

32 33

}

34 35 36 37

public void finalize() { out.close(); }

38 39 40 41 42 43 44 45

public void writeStatus(CarList carList) { out.print(carList.time + " "); for (int i = 0; i < carList.getCount(); i++) { out.print(carList.cars[i].x + " "); } out.println(); } }

7.2

Knihovna Utilities Listing 7.5: conversions.java

1

package utilities;

2 3

public class Conversions {

4 5 6 7

public static double Radians(int degrees) { return degrees*Math.PI / 180; }

8 9

public static double KmHtoMS(double kmh) {

51

return kmh / 3.6;

10

}

11 12

public static double MStoKmH(double ms) { return ms * 3.6; }

13 14 15 16

public static String LaneToStr(int lane) { return (lane==1)?"right lane":"left lane"; }

17 18 19 20

}

Listing 7.6: vectors.java 1

package utilities;

2 3

public class Vectors {

4

public static double[] Vector(double[] A, double[] B) { double[] v = new double[A.length]; for (int i=0; i
5 6 7 8 9 10 11

public static double[] Multiply(double[] v, double k) { double[] newV = new double[v.length]; for (int i = 0; i < v.length; i++) newV[i] = k * v[i]; return newV; }

12 13 14 15 16 17 18

public static double[] LinComb( double alpha, double[] u, double beta, double[] v) { double[] newV = new double[u.length]; for (int i = 0; i < u.length; i++) newV[i] = alpha * u[i] + beta * v[i]; return newV; }

19 20 21 22 23 24 25

public static double[] LinComb( double[] r0, double alpha, double[] u, double beta, double[] v) { double[] newV = new double[u.length]; for (int i = 0; i < u.length; i++) newV[i] = r0[i] + alpha * u[i] + beta * v[i]; return newV; }

26 27 28 29 30 31 32 33

}

7.3

Knihovna mathcomps Listing 7.7: JDrawPanel.java

52

1

package mathcomps;

2 3

import utilities.Vectors;

4 5 6

import javax.swing.*; import java.awt.*;

7 8

public class JDrawPanel extends JPanel {

9 10 11 12 13 14

protected double[] O; protected double x0, y0, x1, y1; protected int w, h; protected double[] u, v; protected boolean initialized = false;

15 16 17 18

public JDrawPanel () { setBackground(Color.BLACK); }

19 20 21 22 23 24 25 26 27 28 29 30 31 32 33

public void setScale (double[] Origin, double[]A, double[]B, double xO, double yO, double xA, double yB) { O = Origin; u = Vectors.Multiply(Vectors.Vector(O, A), 1/(xA-xO) ); v = Vectors.Multiply(Vectors.Vector(O, B), 1/(yB-yO) ); //(10,7) x0 = xO; y0 = yO; x1 = xA; y1 = yB; w = getWidth(); h = getHeight(); initialized = true; } public boolean isInitialized() { return initialized; }

34 35 36 37 38 39

public int[] transform(double x, double y) { double [] vec = Vectors.LinComb(O, x-x0, u, y-y0, v); int [] P = {(int)vec[0], (int)vec[1]}; return P; }

40 41

private int[] lastPoint = new int[2];

42 43 44 45

public void move(Graphics g, double x, double y) { lastPoint = transform(x, y); }

46 47 48 49 50 51 52

public void line(Graphics g, double x, double y) { int[] P = transform(x, y); g.drawLine( lastPoint[0], lastPoint[1], P[0], P[1]); lastPoint[0] = P[0]; lastPoint[1] = P[1]; }

53

53

54 55 56 57 58 59 60 61 62 63 64 65 66 67

public void fillRhomboid( Graphics g, double a1, double b1, double a2, double b2, double a3, double b3, double a4, double b4) { int[] A1 = transform(a1, b1); int[] A2 = transform(a2, b2); int[] A3 = transform(a3, b3); int[] A4 = transform(a4, b4); int[] xs = {A1[0], A2[0], A3[0], A4[0]}; int[] ys = {A1[1], A2[1], A3[1], A4[1]}; Polygon p = new Polygon(xs, ys, 4); g.fillPolygon(p); } }

54

8. Simulace Výstupem programu je tedy soubor, ve kterém je uveden aktuální čas, rychlost a poloha každého vozidla. Tento výstup je klíčový proto, aby se simulace vozidel dala dále použít. Pro zobrazení výsledků se použije naprogramovaný sešit v programu Mathematica. Ten umožňuje zobrazit vozidla jako tok v průběhu času, průběhy drah jednotlivých vozidel, dále pak průběh počtu vozidel v jednotlivých úsecích. Díky tomuto zobrazení se dají snadno vyčíst makroskopické ukazatele a to hlavně hustota vozidel. Parametry všech simulací mají tyto výchozí hodnoty • zrychlení a = 0, 3m/s

2 2

• brzdová konstanta b = 3m/s

• minimální vzdálenost mezi vozidly s0 = 7m • minimální časový odstup mezi vozidly T = 1, 5s • exponent zrychlení δ = 4 Ostatní parametry se budou s ohledem na typ simulace měnit, proto zde nejsou uvedeny.

8.1

Stacionární stav

Stacionární stav je takový stav, kdy jednotlivá vozidla jedou za sebou konstantní rychlostí a mezi sebou udržují konstantní vzdálenosti. K nalezení konstantních rychlostí slouží funkce equilibriumVelocity. Jejími parametry jsou horní a dolní mez rychlosti a vzdálenost mezi vozidly. Tato vzdálenost musí být pro každá dvě vozidla za sebou stejná. Horní mez rychlosti je pak rovna preferované rychlosti vozidla a dolní mez nule. Znamená to tedy, že konstantní rychlost vozidla se bude stanovena mezi těmito mezemi a vozidlo nikdy nedosáhne své preferované rychlosti. Hledá se tedy řešení rovnice " v˙ n = a 1 −



vn vnfree

δ #

 − a

s0 + vn T +

vn (vn −vn−1 ) √ 2 ab

sn

2  = 0,

(8.1)

aby se našel stacionární stav. Po nastavení vfree tedy zůstane pouze jedna neznámá vn . Jediná změna nastává u úplně prvního vozidla, protože před sebou nemá vůbec nikoho a tím pádem je interakční člen roven nule, aint = 0 a rovnice se zjednoduší na "  δ # vn v˙ n = a 1 − = 0, (8.2) vnfree z čehož plyne ze jeho rychlost vn na začátku se musí právě rovnat jeho preferované rychlosti vfree . Pro simulaci stacionárního stavu se použily následující parametry: preferovaná rychlost vozidla vnfree = 120km/h a vzdálenost mezi vozidly sn = 100m. Rychlost, která bude pro všechny vozidla konstantní, potom vychází vn = 110km/h. U cyklické dálnice, kde se poloha vozidla po dosažení konce dálnice vrátí na začátek (jeho poloha po návratu bude nulová), je délka dálnice nastavena na hodnotu 8000m.

55

n@-D 5

n@-D 5

4

4

3

3

2

2

1

1

0

2000

4000

6000

x@mD 10 000

8000

0 4000

6000

a) Počáteční distribuce n@-D 5

4

4

3

3

2

2

1

1

2000

4000

10 000

12 000

x@mD 14 000

b) Konečná distribuce

n@-D 5

0

8000

x@mD 8000

6000

0

2000

c) Počáteční distribuce

4000

6000

x@mD 8000

d) Konečná distribuce

Obrázek 8.1: Simulace stacionárního stavu. a) počáteční distribuce vozidel při nekonečné dálnici b) finální distribuce vozidel při nekonečné dálnici c) počáteční distribuce při cyklické dálnici délky 8 km d) finální distribuce pro cyklickou dálnici délky 8 km

2000

x@mD

1500

1000

500

0

0

100

200

300

400

500

t@sD

Obrázek 8.2: Závislost polohy jednotlivých vozidel na čase v případě stacionárního stavu.

56

12 000 12 000 10 000 8000

8000

x@mD

x@mD

10 000

6000

6000

4000

4000

2000

2000

0

0 0

500

1000

1500

2000

0

t@sD

500

1000

1500

2000

t@sD

a) Nekonečná dálnice

b) Cyklická dálnice

Obrázek 8.3: Závislost poloh jednotlivých vozidel na čase pro celou dobu simulace. Předjíždění v tomto případě musí být zakázáno,1 , aby jednotlivá vozidla neměla potřebu měnit jízdní pruh. Na začátku je tedy hustota rovnoměrně rozložena, a to jak pro nekonečnou dálnici, obrázek 8.1a, tak i pro cyklickou dálnici, obrázek 8.1c. Vozidla jsou od sebe vzdálena konstantně, obrázek 8.2. Z výřezu znázorněném na obrázku 8.2 je také vidět, že použitá numerická metoda v případě stacionárního počátečního rozdělení zachovává ekvidistantní distribuci vozidel. Délky úseků, na kterých je hustota měřena, činí 200 m, proto při vzdálenosti vozidel 100 m je počet vozidel n roven právě dvěma. Konec simulace je znázorněn na obrázcích 8.1b (nekonečná dálnice) a 8.1d (cyklická dálnice). Je vidět, že rozložení automobilů se nemění, pouze se rovnoměrně posouvá rychlostí vn , jež je stejná pro všechna vozidla. Při shlédnutí obrázku 8.3a, kde se na osu x vynáší čas a na ose y figuruje poloha vozidla, je vidět, že se každé jednotlivé vozidlo pohybuje lineárně. Každá úsečka pak odpovídá jednomu vozidlu, které vychází ze své počáteční polohy až do konce simulace. Je tedy vidět, že rovnice (8.1) je opravdu splněna a vozidla se pohybují konstantní rychlostí; zrychlení každého vozidla je nulové. Nic na tom nezmění ani cyklická dálnice, obrázek 8.3b.

8.2

Náhodné rychlosti

Cílem této simulace je ukázat, že při zakázaném předjíždění se vozidla při náhodně nastavených počátečních parametrech dostanou opět do stacionárního stavu. Cílem každého vozidla bude dostat se na svou preferovanou rychlost. Pokud je zakázáno předjíždění, je jasné, že pokud vozidlo bude mít vyšší preferovanou rychlost než vozidlo před ním, na svou preferovanou rychlost by se nemělo dostat a automobil pojede rychlostí nižší, která je rovna preferované rychlosti automobilu před ním. Vozidla při této simulaci budou mít na začátku nastavenou náhodnou preferovanou i počáteční rychlost. Rychlosti budou generovány z Gaussova rozdělení symetrického kolem hodnoty 100km/h. Vygeneruje se číslo z tohoto rozdělení, ke které se poté vynásobí dvaceti a přičte se k němu 100. Může tedy nastat i případ, že auto bude mít vyšší rychlost, než je povolená rychlost na dálnicích v ČR. To ovšem nevadí, protože takoví řidiči se také vyskytují. Pokud by ovšem rychlost vozidla klesla pod 80km/h, pak se rychlost nastaví právě na 80km/h, což odpovídá minimální rychlosti na dálnicích v České republice. Počáteční pozice vozidla se určuje na základě vygenerovaného čísla v rozsahu 0-1. Klasicky zůstávají vzdálenosti mezi vozidly 100m a k tomu se od vygenerovaného čísla odečte 0.5 a 1 Nebo

nastavena taková hodnota parametru p a br , aby pro vozidlo bylo výhodnější jet ve stávajícím pruhu.

57

8000

15 000

6000

x@mD

x@mD

20 000

10 000

5000

4000

2000

0 0

500

1000

1500

2000

2500

3000

0

3500

0

1000

t@sD

2000

3000

4000

5000

t@sD

a) Nekonečná dálnice

b) Cyklická dálnice

Obrázek 8.4: Závislost poloh jednotlivých vozidel na čase při počáteční náhodné rychlosti a poloze pro celou dobu simulace. vynásobí padesáti. Tím pádem se dostane interval +/- 50m. Po prozkoumání obrázku 8.6, který je výřezem obrázku 8.4a je vidět, že startovní pozice opravdu nejsou rovnoměrně rozděleny a vozidla od sebe nejsou konstantně vzdálena. Při pozorování jak se vozidla navzájem dojíždějí začíná být závislost polohy na čase být lineární. Může se stát, že vygenerovaná rychlost je velmi vysoká, zatímco preferovaná může být i 80km/h. Proto je potřeba nechat běžet simulaci potřebně dlouhou dobu, aby se rychlosti jednotlivých vozidel ustálily. To je potom vidět v čase t = 3500s na obrázku 8.4, kde se vozidla shlukla do skupinek, kde první vozidlo brzdí ostatní a ty se tím pádem nemohou dostat na svoji preferovanou rychlost. Jakmile se takto vozidla shlukla, potom platí, že čím déle by se nechala puštěna simulace, tím více by se vozidla přibližovala do stacionárního stavu. Na obrázku 8.4a jsou vidět 3 skupiny vozidel, kde všechny z nich se blíží stacionárnímu stavu. Obrázek pro cyklickou dálnici 8.4b vypadá podobně až na to, že skupinek je podstatně více. To může být způsobeno tím, že vozidlo, které má vygenerovanou nejvzdálenější polohu od počátku (cca. 5000m) interaguje s posledním vozidlem (poloha 0m) nebo jednoduše proto, že byla vygenerovány jiné preferované rychlosti vozidel a jejich polohy. Hustota je na startu simulace přibližně stejná pro nekonečnou dálnici, obrázek 8.5a, tak i pro cyklickou, obrázek 8.5c. Zajímavější obrázky se získávají při zastavení simulace. U nekonečné dálnice je vytvoření skupinek interpretováno zvýšenou hustotou v určitých úsecích, viz obrázek 8.5b. Šířka každého sloupce opět odpovídá hodnotě 200m. U cyklické dálnice o délce 8000m nejsou tyto píky tak patrné, což je způsobeno větším počtem skupinek (obrázek 8.5d) a navíc je dálnice omezena svoji délkou. Simulace ovšem splnila svůj účel, když potvrdila, že vozidla se po určité době dostávají do stacionárního stavu, pokud není povoleno předjíždění.

8.3

Změna parametrů prvního vozidla

V této simulaci půjde o to zkoumat dopravní proud, tz. i stacionární stav systému, při náhlé změně parametrů prvního vozidla. V první fázi se vygenerují vozidla stejně jako v sekci 8.1, obrázek 8.7a. Po nějaké době jízdy se prvnímu vozidlu nastaví jeho preferovaná rychlost na nízkou hodnotu, v této simulaci na 50km/h. Posléze, až se z prvního vozidla stane poslední, se jeho rychlost tentokrát změní na vysokou hodnotu, v tomto případě na 250km/h. Předjíždění musí být povoleno, jeho parametry budou nastaveny standardně, tj. zdvořilostní koeficient p = 0.5

58

n@-D 8

n@-D 8

6

6

4

4

2

2

0

2000

4000

6000

8000

10 000

12 000

x@mD 14 000

0

2000

a) Počáteční distribuce t@sD 5

4

4

3

3

2

2

1

1

0

2000

4000

6000

8000

10 000

12 000

6000

8000

10 000

12 000

x@mD 14 000

12 000

n@-D 14 000

b) Konečná distribuce

t@sD 5

0

4000

n@-D 14 000

0

0

2000

c) Počáteční distribuce

4000

6000

8000

10 000

d) Konečná distribuce

Obrázek 8.5: Simulace při náhodných počátečních polohách a rychlostech. a) počáteční distribuce vozidel při nekonečné dálnici b) finální distribuce při nekonečné dálnici c) počáteční distribuce při cyklické dálnici délky 8 km d) finální distribuce pro cyklickou dálnici délky 8 km

2000

x@mD

1500

1000

500

0

0

200

400

600

800

1000

t@sD

Obrázek 8.6: Náhodné rychlosti - závislost polohy na čase v případě nekonečné dálnice - výřez

59

n@-D 5

n@-D 5

4

4

3

3

2

2

1

1

0

2000

4000

6000

8000

10 000

12 000

x@mD 14 000

a) Počáteční distribuce

0 6000

n@-D 5

4

4

3

3

2

2

1

1

12 000

14 000

16 000

18 000

20 000

22 000

10 000

12 000

14 000

16 000

18 000

x@mD 20 000

b) Distribuce, rychlost prvního vozidla 50km/h

n@-D 5

0 10 000

8000

24 000

x@mD

0 20 000

c) Distribuce po předjíždění

25 000

30 000

35 000

40 000

x@mD 45 000

d) Konečná distribuce

Obrázek 8.7: Simulace při snížení a následném zvýšení rychlosti čelního vozidla. a) počáteční distribuce vozidel při nekonečné dálnici b) distribuce při rychlosti čelního vozidla 50 km/h c) distribuce poté po předjetí čelního vozidla d) finální distribuce a preference pravého jízdního pruhu. Předpoklad je takový, že vozidla začnou plynule vozidlo předjíždět. Poté, co se mu nastaví velmi vysoká rychlost, bude zase vozidlo plynule předjíždět všechna vozidla a systém až na toto vozidlo, které mění svojí preferovanou rychlost, zůstanou ve stacionárním stavu. Jak je vidět u nekonečně dálnice, výsledky jsou poměrně překvapující. Při snížení rychlosti prvního vozidla se ještě nic neděje, automobil začne klasicky zpomalovat na svou preferovanou rychlost 50km/h a ostatní vozidla ho začnou klasicky předjíždět. Ovšem v tu chvíli už se začíná pomalu nabourávat stabilita systému. Je to způsobeno tím, že vozidlo začne předjíždět až ve chvíli, kdy mu to umožní podmínka motivačního kritéria, viz odstavec 5.2.1. Jeho rychlost se mezitím sníží a on až poté začne předjíždět. Tento proces se šíří mezi jednotlivými vozidly. Poté co ho všechna předjedou, tak teď už poslední vozidlo zrychlí na preferovanou rychlost 250km/h. Vozidla před ním ovšem nejsou ve stacionárním stavu, došlo tam k narušení stability a vozidla před ním se navzájem předjíždějí, takže i když má vozidlo nastavenou preferovanou rychlost velmi vysokou, není šance aby se před tyto vozidla dostal. Znázorněno je to na obrázku 8.7d, kde vozidlo předjelo asi jen 6 vozidel oproti tomu, když ho všechna vozidla předjela, a to za stejný čas. Jeho konečná poloha se nachází v růžově zvýrazněném sloupci vozidel. Tato simulace ukazuje vysokou závislost na nastavení parametrů nejen předjížděcího modelu, ale také na nastavení chvíle, kdy řidič začne opět zrychlovat. Na obrázku 8.9 je potom vidět závislost dráhy na čase. Je vidět, že čelní vozidlo je postupně předjížděno ostatními, na své původní místo se ovšem za stejný čas nevrátí a skončí někde mezi

60

22 000

30 000

21 000 29 000

19 000

x@mD

x@mD

20 000

18 000

28 000

27 000

17 000 26 000 16 000 15 000 5000

5500

6000

6500

25 000 5000

7000

5500

t@sD

6000

6500

7000

t@sD

a) výřez - předjíždění posledního vozidla

b) výřez - posledních 2000 s

Obrázek 8.8: Simulace při snížení a následném zvýšení rychlosti čelního vozidla. a) závislost polohy vozidel na čase - zvýšení rychlosti posledního vozidla b) závislost polohy vozidel na čase - posledních 2000 s

40 000

x@mD

30 000

20 000

10 000

0 0

1000

2000

3000

4000

5000

6000

7000

t@sD Obrázek 8.9: Simulace při snížení a následném zvýšení rychlosti čelního vozidla - závislost polohy na čase - nekonečná dálnice

61

8000

x@mD

6000

4000

2000

0

0

1000

2000

3000

4000

5000

t@sD Obrázek 8.10: Simulace při snížení a následném zvýšení rychlosti čelního vozidla - závislost polohy na čase - cyklická dálnice ostatními vozidly, což dokumentuje obrázek 8.8a. Na výřezu 8.8b je potom vidět, že všechna vozidla nejsou ve stacionárním stavu, že čelní vozidlo celkově stabilitu systému narušilo. U cyklické dálnice se stacionární stav nenaboural tolik, jako u nekonečné, je pouze vidět zvětšení vzdálenosti mezi jednotlivými vozidly při zpomalení vozidla, viz obrázek 8.10. Není úplně přesně jasné, proč se dálnice v případě cykličnosti chová jinak než u nekonečné. Hustota u cyklické dálnice odpovídá počátečním předpokladům, vozidlo předjíždí plynule ostatní vozidla ve stacionárním stavu, viz obrázek 8.11. Růžově je zobrazen sloupec s vozidly, ve kterém se nachází předjížděcí.

8.4

Uzavření pruhu komunikace

V poslední simulaci se v určitém úseku uzavře pravý jízdní pruh komunikace. Vozidla v tomto případě budou mít na začátku nastavenu náhodnou rychlost a polohu stejně jako v simulaci 8.2. To znamená, že pokud by zde uzavírka nebyla, vozidla postupně přejdou do skupin stacionárního stavu. Hlavní problém této simulace se týká nastavení parametrů předjížděcího modelu. Nejprve byla vyzkoušena simulace s parametry používanými do teď, to jest zdvořilostní koeficient p = 0.5, dále br = 0.3 pro vozidla v pravém jízdním pruhu a br = −0.1 pro levý pruh, což značí, že řidič se ohlíží na vozidlo jedoucí ve vedlejším pruhu, pokud je v levém pruhu, snaží se hned zařadit do pravého. Takto spuštěná simulace vůbec neodpovídala realitě. Vozidla, která měla nízkou preferovanou rychlost a jela v pravém jízdním pruhu vůbec neměla šanci se této uzavírce pruhu vyhnout. Je to způsobeno tím, že berou přílišný ohled na vozidlo jedoucí v levém pruhu. Proto bylo nutné tyto parametry přenastavit. Čím blíže se bude vozidlo blížit uzavírce, tím méně ho zajímá vozidlo ve vedlejším pruhu. Přeneseno do předjížděcího modelu, 500 m před překážkou má každé vozidlo nastavenou hodnotu zdvořilostního koeficientu na hodnotu p = 0.2 a 150 m už jen p = 0.05. Potom, co se uzavírka pruhu objede, se opět nastaví p na výchozí hodnotu 0.5. Preference levého či pravého pruhu je

62

n@-D 6

5

4

3

2

1

0

2000

4000

6000

8000

x@mD

Obrázek 8.11: Simulace při snížení a následném zvýšení rychlosti čelního vozidla - finální distribuce - cyklická dálnice pak nastavena na nulu, br = 0. Vozidla byla vygenerována až do polohy 5000m a jízdní pruh byl uzavřen od 6050m do 6150m (na obrázku 8.12a znázorněno dvěma vozidly v úseku 6000 - 6200m). Za tímto místem už by se vozidla měla snažit klasicky dostat na svoji preferovanou rychlost. Na obrázku 8.12a je vidět rozložení hustoty vozidel po částečně náhodném vygenerování jejich pozic. Pokud simulace poběží dostatečně dlouho (více než tři hodiny) výsledkem je obrázek 8.12b, kde jsou znatelně vidět jednotlivé skupiny vozidel, které se shlukují. Zaprvé je to tím, že se k sobě shlukují vozidla se stejnými preferovanými rychlostmi a za druhé, že po předjetí překážky vozidla již nejeví takový zájem předjíždět vozidla před sebou, i když mají třeba o několik km nižší preferovanou rychlost. Průběh polohy na čase je potom znázorněn na obrázku 8.13a. Je vidět, že u překážky na 6000m (na ose y) vozidla velmi zpomalují, navíc některé musejí čekat, než přejede vozidlo za ním ve druhém jízdním pruhu a až pak se mohou zařadit a přejet si. Druhý extrém nastává v případě, že vozidlo v levém jízdním pruhu jede rychle, a vozidlo mu tam zprava „vjedeÿ. Nastává rychlé brždění, opět znázorněno na obrázku 8.13a. U cyklické dálnice se vygenerovala opět vozidla na prvních 5000m, viz. obrázek 8.12c . Výsledky simulace docela překvapují, protože vozidla opět utváří určité skupinky, viz. obrázek 8.12d, což by se na tak krátké vzdálenosti ani nepředpokládalo. Na obrázku 8.13b je znázorněna závislost polohy na čase. Některá vozidla pouze musejí přibrzdit. Poté, co se však utvoří skupinky, již není nutné při předjíždění překážky zpomalovat.

63

n@-D 5

n@-D 12

10

4

8 3 6 2 4 1

2

0

2000

4000

6000

8000

10 000

x@mD 14 000

12 000

0 20 000

25 000

a) Počáteční distribuce n@-D 6

5

5

4

4

3

3

2

2

1

1

2000

4000

6000

x@mD 40 000

35 000

b) Konečná distribuce

n@-D 6

0

30 000

8000

x@mD

0

2000

c) Počáteční distribuce

4000

6000

8000

x@mD

d) Konečná distribuce

Obrázek 8.12: Simulace uzavření jízdního pruhu. a) počáteční distribuce vozidel při nekonečné dálnici b) finální distribuce při nekonečné dálnici c) počáteční distribuce při cyklické dálnici délky 8 km d) finální distribuce pro cyklickou dálnici při cyklické dálnici délky 8 km

60 000 8000 50 000 6000

x@mD

x@mD

40 000 30 000

4000

20 000 2000 10 000 0 0

2000

4000

6000

8000

0

10 000

t@sD

0

2000

4000

6000

8000

10 000

12 000

t@sD

a) Nekonečná dálnice

b) Cyklická dálnice

Obrázek 8.13: Simulace při uzavření jednoho pruhu a) závislost polohy vozidel na čase - nekonečná dálnice b) závislost polohy vozidel na čase - cyklická dálnice

64

Závěr Jak už je zmíněno v úvodu, cílem této práce nebylo dosažení nových vědeckých poznatků, nýbrž seznámení se s problematikou dopravních modelů a jejich programování. Za tímto účelem byly popsány dva základní typy modelů: makroskopické a mikroskopické. Z vlastností těchto typů modelů vyplynulo, že pro další studium bude vhodné se primárně zabývat mikroskopickými modely. Tyto modely jsou převážně popsány diferenciálními rovnicemi. Porozumění jejich tematice je podmíněno znalostmi z oblasti teorie stability. Tímto se zabývá čtvrtá kapitola této práce. Pro realizaci dopravních modelů založených na diferenciálních rovnicích je nutno použít numerické postupy a algoritmy. Studovali jsme základní Eulerovu metodu numerického řešení diferenciálních rovnic. Poté jsme implementovali metody vyšších řádu, konkrétně metody Runge-Kutta druhého a čtvrtého řádu. Ukázalo se však, že všechny tyto metody poskytují pro studovaný problém identické výsledky. Na základě teoretického popisu modelu IDM byla vytvořena základní kostra programu. V této podobě program pracoval pouze s číselnými vstupy a výstupy. Díky tomu mohl být vyzkoušen matematický popis tohoto modelu a mohla být ověřena jeho funkčnost. Po úspěšných výsledcích testování bylo k programu přidáno grafické prostředí, v němž bylo možné simulovat i složitější problémy. Pro implementaci byl použit programovací jazyk Java. Po implemntaci IDM byla aplikace postupně rozšiřována o předjížděcí model MOBIL. Naprogramování tohoto modelu bylo jednoduché, složitější však bylo zjištění jeho správných vstupních parametrů. Dosažení adekvátních výsledků vyžadovalo dlouhé testování spojené s mnohonásobným přenastavováním jednotlivých hodnot. Po těchto úspěšně zakončených pokusech se přistoupilo k samotným simulacím. Pro simulaci stacionárního stavu bylo třeba naprogramovat funkci, pomocí které by bylo možno spočítat jednotlivé kritické body a a díky nim stanovit popis tohoto stavu. K jejímu vytvoření byla použita metoda sečen. Následně byly provedeny simulace s využitím modelu IDM (a rozšířeného modelu IDM+MOBIL) pro tento stav a pro další dopravní simulace. Výsledky simulací bylo výhodné popsat dvojím způsobem: makroskopicky v pojmech distribucí počtu vozidel (mj. hustoty) a mikroskopicky se zaměřením na jednotlivá vozidla. Nejprve byla otestována zvolená numerická metoda na stacionárním stavu modelu IDM, kde výsledky simulace potvrdily vhodnost jejího použití. Postupně byly testovány tři případy na nekonečné i cyklické dálnici. Zjistili jsme, že ve všech případech mají vozidla tendenci dostávat se na svoji preferovanou rychlost, ale s ohledem na ostatní vozidla tendují ke stacionárnímu stavu, přičemž u třetího případu (podkapitola 8.3) by stacionárního stavu bylo dosaženo za poměrně dlouhý časový úsek. Výzkum v této oblasti je možno rozšiřovat směrem k dalším mikroskopickým modelům či zvolit makroskopický model pro hustší provoz na pozemní komunikaci.

65

Literatura [1] Barceló J, Fundamentals of traffic simulation, proceedings, International Series in Operations Research and management science, Springer, 2010 [2] Barrow J, Nové teorie všeho, Dokořán, Praha, 2008 [3] Brdička M, Samek L, Sopko B, Mechanika kontinua, Academia, Praha, 2011 [4] Haberman R, Mathematical models: Mechanical vibrations, population dynamics and traffic flow, Society for Industrial and Applied Mathematics, Philadelphia, 1988 [5] Helbing D, Herrmann H J, Schreckenberg M, Wolf D E, Microscopic Simulation of Congested Traffic v knížce Traffic and Granular Flow, Springer, Berlin, 2000 [6] Horák J, Krlín L, Raidl A, Deterministický chaos a jeho fyzikální aplikace, Academia 2003 [7] Xiaoliang Ma, A Neural-Fuzzy Framework for Modeling Car-following Behavior [online], dostupné z http://www.ctr.kth.se/publications/ctr2006_08.pdf [8] May A D, Traffic flow fundamentals, Prentice Hall, 1989 [9] Přikryl P, Numerické metody matematické analýzy, SNTL, Praha, 1988 [10] Scholtz M, Classical mechanics and deterministic chaos [online], dostupné z http://www. fd.cvut.cz/personal/scholma1/ [11] Scholtz M, Vaniš M, Veselý P, Matějka P, Applied mathematics on Faculty of Transportation Sciences, vyjde ve sborníku k výročí Fakulty dopravní ČVUT [12] Treiber M, Hennecke A, Helbing D, Congested traffic states in empirical observations and microscopic simulations, Physical Review E, 62 (2), pp. 1805–1824, 2000 [13] Treiber M, Microsimulation of road traffic flow [online], dostupné z http://www. traffic-simulation.de/ [14] Vitásek E, Numerické metody, SNTL, Praha, 1987

66

Seznam obrázků Fundamentální diagramy pro Greenshieldsův model. a) Závislost v = v(k). b) Závislost q = q(k). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

2.1

Zóny pro reakci řidičů v závislosti na rozdílu rychlostí a poloh . . . . . . . . . . .

13

3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8

Lorenzův atraktor. . . . . . . Trajektorie pro λ1,2 > 0 . . . Trajektorie pro λ1,2 < 0 . . . Trajektorie pro λ1 · λ2 < 0 . . Trajektorie pro α = 0, β < 0 . Trajektorie pro α = 0, β > 0 . Trajektorie pro α < 0, β > 0 . Trajektorie pro α > 0, β > 0 .

. . . . . . . .

15 16 16 17 17 17 18 18

4.1

Grafické znázornění numerických metod . . . . . . . . . . . . . . . . . . . . . . .

21

5.1 5.2

Proměnné jednotlivých aut . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Schéma situace při výpočtu předjíždění . . . . . . . . . . . . . . . . . . . . . . .

25 26

6.1 6.2 6.3

Snímek obrazovky programu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Popis bodu pomocí lineární kombinace . . . . . . . . . . . . . . . . . . . . . . . . Konstanty použité pro vykreslení dálnice . . . . . . . . . . . . . . . . . . . . . . .

30 31 36

1.1

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

8.1

Simulace stacionárního stavu. a) počáteční distribuce vozidel při nekonečné dálnici b) finální distribuce vozidel při nekonečné dálnici c) počáteční distribuce při cyklické dálnici délky 8 km d) finální distribuce pro cyklickou dálnici délky 8 km 8.2 Závislost polohy jednotlivých vozidel na čase v případě stacionárního stavu. . . . 8.3 Závislost poloh jednotlivých vozidel na čase pro celou dobu simulace. . . . . . . . 8.4 Závislost poloh jednotlivých vozidel na čase při počáteční náhodné rychlosti a poloze pro celou dobu simulace. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.5 Simulace při náhodných počátečních polohách a rychlostech. a) počáteční distribuce vozidel při nekonečné dálnici b) finální distribuce při nekonečné dálnici c) počáteční distribuce při cyklické dálnici délky 8 km d) finální distribuce pro cyklickou dálnici délky 8 km . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.6 Náhodné rychlosti - závislost polohy na čase v případě nekonečné dálnice - výřez 8.7 Simulace při snížení a následném zvýšení rychlosti čelního vozidla. a) počáteční distribuce vozidel při nekonečné dálnici b) distribuce při rychlosti čelního vozidla 50 km/h c) distribuce poté po předjetí čelního vozidla d) finální distribuce . . . 8.8 Simulace při snížení a následném zvýšení rychlosti čelního vozidla. a) závislost polohy vozidel na čase - zvýšení rychlosti posledního vozidla b) závislost polohy vozidel na čase - posledních 2000 s . . . . . . . . . . . . . . . . . . . . . . . . . . 8.9 Simulace při snížení a následném zvýšení rychlosti čelního vozidla - závislost polohy na čase - nekonečná dálnice . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.10 Simulace při snížení a následném zvýšení rychlosti čelního vozidla - závislost polohy na čase - cyklická dálnice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.11 Simulace při snížení a následném zvýšení rychlosti čelního vozidla - finální distribuce - cyklická dálnice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

56 56 57 58

59 59

60

61 61 62 63

8.12 Simulace uzavření jízdního pruhu. a) počáteční distribuce vozidel při nekonečné dálnici b) finální distribuce při nekonečné dálnici c) počáteční distribuce při cyklické dálnici délky 8 km d) finální distribuce pro cyklickou dálnici při cyklické dálnici délky 8 km . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.13 Simulace při uzavření jednoho pruhu a) závislost polohy vozidel na čase - nekonečná dálnice b) závislost polohy vozidel na čase - cyklická dálnice . . . . . . . .

68

64 64

Seznam tabulek 5.1

Hodnoty používané pro simulaci Zdroj:[13] . . . . . . . . . . . . . . . . . . . . . .

69

25