Line data Source code
1 : /****************************************************************************/
2 : // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
3 : // Copyright (C) 2001-2026 German Aerospace Center (DLR) and others.
4 : // This program and the accompanying materials are made available under the
5 : // terms of the Eclipse Public License 2.0 which is available at
6 : // https://www.eclipse.org/legal/epl-2.0/
7 : // This Source Code may also be made available under the following Secondary
8 : // Licenses when the conditions for such availability set forth in the Eclipse
9 : // Public License 2.0 are satisfied: GNU General Public License, version 2
10 : // or later which is available at
11 : // https://www.gnu.org/licenses/old-licenses/gpl-2.0-standalone.html
12 : // SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
13 : /****************************************************************************/
14 : /// @file ROMAAssignments.cpp
15 : /// @author Yun-Pang Floetteroed
16 : /// @author Laura Bieker
17 : /// @author Michael Behrisch
18 : /// @date Feb 2013
19 : ///
20 : // Assignment methods
21 : /****************************************************************************/
22 : #include <config.h>
23 :
24 : #include <vector>
25 : #include <algorithm>
26 : #include <utils/common/SUMOTime.h>
27 : #include <utils/distribution/Distribution_Points.h>
28 : #include <utils/router/RouteCostCalculator.h>
29 : #include <utils/router/SUMOAbstractRouter.h>
30 : #include <router/ROEdge.h>
31 : #include <router/RONet.h>
32 : #include <router/RORoute.h>
33 : #include <od/ODMatrix.h>
34 : #include "ROMAEdge.h"
35 : #include "ROMAAssignments.h"
36 :
37 :
38 : // ===========================================================================
39 : // static member variables
40 : // ===========================================================================
41 : std::map<const ROEdge* const, double> ROMAAssignments::myPenalties;
42 :
43 :
44 : // ===========================================================================
45 : // method definitions
46 : // ===========================================================================
47 :
48 73 : ROMAAssignments::ROMAAssignments(const SUMOTime begin, const SUMOTime end, const bool additiveTraffic,
49 : const double adaptionFactor, const int maxAlternatives, const bool defaultCapacities,
50 : RONet& net, ODMatrix& matrix, SUMOAbstractRouter<ROEdge, ROVehicle>& router,
51 73 : OutputDevice* netloadOutput)
52 73 : : myBegin(begin), myEnd(end), myAdditiveTraffic(additiveTraffic), myAdaptionFactor(adaptionFactor),
53 73 : myMaxAlternatives(maxAlternatives), myUseDefaultCapacities(defaultCapacities),
54 73 : myNet(net), myMatrix(matrix), myRouter(router), myNetloadOutput(netloadOutput) {
55 73 : myDefaultVehicle = new ROVehicle(SUMOVehicleParameter(), nullptr, net.getVehicleTypeSecure(DEFAULT_VTYPE_ID), &net);
56 73 : }
57 :
58 :
59 73 : ROMAAssignments::~ROMAAssignments() {
60 73 : delete myDefaultVehicle;
61 73 : }
62 :
63 : // based on the definitions in PTV-Validate and in the VISUM-Cologne network
64 : double
65 162597 : ROMAAssignments::getCapacity(const ROEdge* edge) const {
66 162597 : if (edge->isTazConnector()) {
67 : return 0;
68 : }
69 162597 : const int roadClass = myUseDefaultCapacities ? -1 : -edge->getPriority();
70 : // TODO: differ road class 1 from the unknown road class 1!!!
71 162597 : if (edge->getNumLanes() == 0) {
72 : // TAZ have no cost
73 : return 0;
74 162597 : } else if (roadClass == 0 || roadClass == 1) {
75 160317 : return edge->getNumLanes() * 2000.; //CR13 in table.py
76 2280 : } else if (roadClass == 2 && edge->getSpeedLimit() <= 11.) {
77 0 : return edge->getNumLanes() * 1333.33; //CR5 in table.py
78 2280 : } else if (roadClass == 2 && edge->getSpeedLimit() > 11. && edge->getSpeedLimit() <= 16.) {
79 0 : return edge->getNumLanes() * 1500.; //CR3 in table.py
80 2280 : } else if (roadClass == 2 && edge->getSpeedLimit() > 16.) {
81 0 : return edge->getNumLanes() * 2000.; //CR13 in table.py
82 2280 : } else if (roadClass == 3 && edge->getSpeedLimit() <= 11.) {
83 0 : return edge->getNumLanes() * 800.; //CR5 in table.py
84 2280 : } else if (roadClass == 3 && edge->getSpeedLimit() > 11. && edge->getSpeedLimit() <= 13.) {
85 0 : return edge->getNumLanes() * 875.; //CR5 in table.py
86 2280 : } else if (roadClass == 3 && edge->getSpeedLimit() > 13. && edge->getSpeedLimit() <= 16.) {
87 0 : return edge->getNumLanes() * 1500.; //CR4 in table.py
88 2280 : } else if (roadClass == 3 && edge->getSpeedLimit() > 16.) {
89 0 : return edge->getNumLanes() * 1800.; //CR13 in table.py
90 2280 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() <= 5.) {
91 0 : return edge->getNumLanes() * 200.; //CR7 in table.py
92 2280 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 5. && edge->getSpeedLimit() <= 7.) {
93 0 : return edge->getNumLanes() * 412.5; //CR7 in table.py
94 2280 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 7. && edge->getSpeedLimit() <= 9.) {
95 0 : return edge->getNumLanes() * 600.; //CR6 in table.py
96 2280 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 9. && edge->getSpeedLimit() <= 11.) {
97 0 : return edge->getNumLanes() * 800.; //CR5 in table.py
98 2280 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 11. && edge->getSpeedLimit() <= 13.) {
99 0 : return edge->getNumLanes() * 1125.; //CR5 in table.py
100 2280 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 13. && edge->getSpeedLimit() <= 16.) {
101 2250 : return edge->getNumLanes() * 1583.; //CR4 in table.py
102 30 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 16. && edge->getSpeedLimit() <= 18.) {
103 0 : return edge->getNumLanes() * 1100.; //CR3 in table.py
104 30 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 18. && edge->getSpeedLimit() <= 22.) {
105 0 : return edge->getNumLanes() * 1200.; //CR3 in table.py
106 30 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 22. && edge->getSpeedLimit() <= 26.) {
107 0 : return edge->getNumLanes() * 1300.; //CR3 in table.py
108 30 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 26.) {
109 0 : return edge->getNumLanes() * 1400.; //CR3 in table.py
110 : }
111 30 : return edge->getNumLanes() * 800.; //CR5 in table.py
112 : }
113 :
114 :
115 : // based on the definitions in PTV-Validate and in the VISUM-Cologne network
116 : double
117 242600 : ROMAAssignments::capacityConstraintFunction(const ROEdge* edge, const double flow) const {
118 242600 : if (edge->isTazConnector()) {
119 : return 0;
120 : }
121 161378 : const int roadClass = myUseDefaultCapacities ? -1 : -edge->getPriority();
122 161378 : const double capacity = getCapacity(edge);
123 : // TODO: differ road class 1 from the unknown road class 1!!!
124 161378 : if (edge->getNumLanes() == 0) {
125 : // TAZ have no cost
126 : return 0;
127 161378 : } else if (roadClass == 0 || roadClass == 1) {
128 159116 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 1.3)) * 2.); //CR13 in table.py
129 2262 : } else if (roadClass == 2 && edge->getSpeedLimit() <= 11.) {
130 0 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.9)) * 3.); //CR5 in table.py
131 2262 : } else if (roadClass == 2 && edge->getSpeedLimit() > 11. && edge->getSpeedLimit() <= 16.) {
132 0 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 1.)) * 2.); //CR3 in table.py
133 2262 : } else if (roadClass == 2 && edge->getSpeedLimit() > 16.) {
134 0 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 1.3)) * 2.); //CR13 in table.py
135 2262 : } else if (roadClass == 3 && edge->getSpeedLimit() <= 11.) {
136 0 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.9)) * 3.); //CR5 in table.py
137 2262 : } else if (roadClass == 3 && edge->getSpeedLimit() > 11. && edge->getSpeedLimit() <= 13.) {
138 0 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.9)) * 3.); //CR5 in table.py
139 2262 : } else if (roadClass == 3 && edge->getSpeedLimit() > 13. && edge->getSpeedLimit() <= 16.) {
140 0 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.7 * (flow / (capacity * 1.)) * 2.); //CR4 in table.py
141 2262 : } else if (roadClass == 3 && edge->getSpeedLimit() > 16.) {
142 0 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 1.3)) * 2.); //CR13 in table.py
143 2262 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() <= 5.) {
144 0 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.5)) * 3.); //CR7 in table.py
145 2262 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 5. && edge->getSpeedLimit() <= 7.) {
146 0 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.5)) * 3.); //CR7 in table.py
147 2262 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 7. && edge->getSpeedLimit() <= 9.) {
148 0 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.8)) * 3.); //CR6 in table.py
149 2262 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 9. && edge->getSpeedLimit() <= 11.) {
150 0 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.9)) * 3.); //CR5 in table.py
151 2262 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 11. && edge->getSpeedLimit() <= 13.) {
152 0 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.9)) * 3.); //CR5 in table.py
153 2262 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 13. && edge->getSpeedLimit() <= 16.) {
154 2235 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.7 * (flow / (capacity * 1.)) * 2.); //CR4 in table.py
155 27 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 16. && edge->getSpeedLimit() <= 18.) {
156 0 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 1.)) * 2.); //CR3 in table.py
157 27 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 18. && edge->getSpeedLimit() <= 22.) {
158 0 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 1.)) * 2.); //CR3 in table.py
159 27 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 22. && edge->getSpeedLimit() <= 26.) {
160 0 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 1.)) * 2.); //CR3 in table.py
161 27 : } else if ((roadClass >= 4 || roadClass == -1) && edge->getSpeedLimit() > 26.) {
162 0 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 1.)) * 2.); //CR3 in table.py
163 : }
164 27 : return edge->getLength() / edge->getSpeedLimit() * (1. + 1.*(flow / (capacity * 0.9)) * 3.); //CR5 in table.py
165 : }
166 :
167 :
168 : bool
169 40837 : ROMAAssignments::addRoute(const ConstROEdgeVector& edges, std::vector<RORoute*>& paths, std::string routeId, double prob) {
170 : std::vector<RORoute*>::iterator p;
171 48351 : for (p = paths.begin(); p != paths.end(); p++) {
172 45704 : if (edges == (*p)->getEdgeVector()) {
173 : break;
174 : }
175 : }
176 40837 : if (p == paths.end()) {
177 2647 : paths.push_back(new RORoute(routeId, 0., prob, edges, nullptr, StopParVector()));
178 2647 : return true;
179 : }
180 38190 : (*p)->addProbability(prob);
181 : std::iter_swap(paths.end() - 1, p);
182 38190 : return false;
183 : }
184 :
185 :
186 : const ConstROEdgeVector
187 40875 : ROMAAssignments::computePath(ODCell* cell, const SUMOTime time, const double probability, SUMOAbstractRouter<ROEdge, ROVehicle>* router, bool setBulkMode) {
188 81490 : const ROEdge* const from = myNet.getEdge(cell->origin + (cell->originIsEdge ? "" : "-source"));
189 40875 : if (from == nullptr) {
190 0 : throw ProcessError(TLF("Unknown origin '%'.", cell->origin));
191 : }
192 81490 : const ROEdge* const to = myNet.getEdge(cell->destination + (cell->destinationIsEdge ? "" : "-sink"));
193 40875 : if (to == nullptr) {
194 0 : throw ProcessError(TLF("Unknown destination '%'.", cell->destination));
195 : }
196 40875 : ConstROEdgeVector edges;
197 40875 : if (router == nullptr) {
198 15 : router = &myRouter;
199 : }
200 40875 : if (myMaxAlternatives > 0 && (int)cell->pathsVector.size() < myMaxAlternatives) {
201 40837 : router->compute(from, to, myDefaultVehicle, time, edges);
202 40837 : if (setBulkMode) {
203 40762 : router->setBulkMode(true);
204 : }
205 122511 : if (addRoute(edges, cell->pathsVector, cell->origin + cell->destination + toString(cell->pathsVector.size()), probability)) {
206 : return edges;
207 : }
208 : } else {
209 : double minCost = std::numeric_limits<double>::max();
210 : RORoute* minRoute = nullptr;
211 76 : for (RORoute* const p : cell->pathsVector) {
212 38 : const double cost = router->recomputeCosts(edges, myDefaultVehicle, time);
213 38 : if (cost < minCost) {
214 : minCost = cost;
215 : minRoute = p;
216 : }
217 : }
218 38 : minRoute->addProbability(probability);
219 : }
220 38228 : return ConstROEdgeVector();
221 40875 : }
222 :
223 :
224 : void
225 4 : ROMAAssignments::getKPaths(const int kPaths, const double penalty) {
226 8 : for (ODCell* const c : myMatrix.getCells()) {
227 : myPenalties.clear();
228 14 : for (int k = 0; k < kPaths; k++) {
229 67 : for (const ROEdge* const e : computePath(c)) {
230 57 : myPenalties[e] += penalty;
231 10 : }
232 : }
233 : }
234 : myPenalties.clear();
235 4 : }
236 :
237 :
238 : void
239 73 : ROMAAssignments::resetFlows() {
240 73 : const double begin = STEPS2TIME(MIN2(myBegin, myMatrix.getCells().front()->begin));
241 2206 : for (ROEdge* e : ROEdge::getAllEdges()) {
242 : ROMAEdge* const edge = static_cast<ROMAEdge*>(e);
243 2133 : edge->setFlow(begin, STEPS2TIME(myEnd), 0.);
244 2133 : edge->setHelpFlow(begin, STEPS2TIME(myEnd), 0.);
245 : }
246 73 : }
247 :
248 :
249 : void
250 214 : ROMAAssignments::writeInterval(const SUMOTime begin, const SUMOTime end) {
251 214 : if (myNetloadOutput != nullptr) {
252 64 : myNetloadOutput->openTag(SUMO_TAG_INTERVAL).writeAttr(SUMO_ATTR_BEGIN, time2string(begin)).writeAttr(SUMO_ATTR_END, time2string(end));
253 6491 : for (ROEdge* e : ROEdge::getAllEdges()) {
254 : const ROMAEdge* const edge = static_cast<ROMAEdge*>(e);
255 6459 : if (edge->getFunction() == SumoXMLEdgeFunc::NORMAL) {
256 1219 : myNetloadOutput->openTag(SUMO_TAG_EDGE).writeAttr(SUMO_ATTR_ID, edge->getID());
257 1219 : const double traveltime = edge->getTravelTime(getDefaultVehicle(), STEPS2TIME(begin));
258 1219 : const double speed = edge->getLength() / traveltime;
259 1219 : const double flow = edge->getFlow(STEPS2TIME(begin));
260 1219 : const double timeGap = STEPS2TIME(end - begin) / flow;
261 1219 : const double spaceGap = timeGap * speed;
262 1219 : const double density = 1000.0 / spaceGap;
263 1219 : const double laneDensity = density / edge->getNumLanes();
264 1219 : myNetloadOutput->writeAttr(SUMO_ATTR_TRAVELTIME, traveltime);
265 1219 : myNetloadOutput->writeAttr(SUMO_ATTR_SPEED, speed);
266 1219 : myNetloadOutput->writeAttr(SUMO_ATTR_SPEEDREL, speed / edge->getSpeedLimit());
267 1219 : myNetloadOutput->writeAttr(SUMO_ATTR_ENTERED, flow);
268 1219 : myNetloadOutput->writeAttr(SUMO_ATTR_DENSITY, density);
269 1219 : myNetloadOutput->writeAttr(SUMO_ATTR_LANEDENSITY, laneDensity);
270 1219 : myNetloadOutput->writeAttr("flowCapacityRatio", 100. * flow / getCapacity(edge));
271 2438 : myNetloadOutput->closeTag();
272 : }
273 : }
274 64 : myNetloadOutput->closeTag();
275 : }
276 214 : }
277 :
278 :
279 : void
280 69 : ROMAAssignments::incremental(const int numIter, const bool verbose) {
281 : SUMOTime lastBegin = -1;
282 69 : std::vector<int> intervals;
283 69 : int count = 0;
284 2109 : for (const ODCell* const c : myMatrix.getCells()) {
285 2040 : if (c->begin != lastBegin) {
286 210 : intervals.push_back(count);
287 210 : lastBegin = c->begin;
288 : }
289 2040 : count++;
290 : }
291 : lastBegin = -1;
292 279 : for (std::vector<int>::const_iterator offset = intervals.begin(); offset != intervals.end(); offset++) {
293 210 : std::vector<ODCell*>::const_iterator firstCell = myMatrix.getCells().begin() + (*offset);
294 : std::vector<ODCell*>::const_iterator lastCell = myMatrix.getCells().end();
295 210 : if (offset != intervals.end() - 1) {
296 141 : lastCell = myMatrix.getCells().begin() + (*(offset + 1));
297 : }
298 210 : const SUMOTime intervalStart = (*firstCell)->begin;
299 210 : if (verbose) {
300 0 : WRITE_MESSAGE(" starting interval " + time2string(intervalStart));
301 : }
302 : std::map<const ROMAEdge*, double> loadedTravelTimes;
303 27021 : for (ROEdge* e : ROEdge::getAllEdges()) {
304 26811 : const ROMAEdge* const edge = static_cast<ROMAEdge*>(e);
305 26811 : if (edge->hasLoadedTravelTime(STEPS2TIME(intervalStart))) {
306 2024 : loadedTravelTimes[edge] = edge->getTravelTime(myDefaultVehicle, STEPS2TIME(intervalStart));
307 : }
308 : }
309 210 : if (loadedTravelTimes.empty()) {
310 : // we don't have loaded travel times for this interval but we may have set ourselves some for the interval before but no need to warn
311 : ROEdge::disableTimelineWarning();
312 : }
313 4440 : for (int t = 0; t < numIter; t++) {
314 4230 : if (verbose) {
315 0 : WRITE_MESSAGE(" starting iteration " + toString(t));
316 : }
317 4230 : std::string lastOrigin = "";
318 : int workerIndex = 0;
319 45090 : for (std::vector<ODCell*>::const_iterator i = firstCell; i != lastCell; i++) {
320 40860 : ODCell* const c = *i;
321 40860 : const double linkFlow = c->vehicleNumber / numIter;
322 40860 : const SUMOTime begin = myAdditiveTraffic ? myBegin : c->begin;
323 : #ifdef HAVE_FOX
324 40860 : if (myNet.getThreadPool().size() > 0) {
325 60 : if (lastOrigin != c->origin) {
326 60 : workerIndex++;
327 60 : if (workerIndex == myNet.getThreadPool().size()) {
328 : workerIndex = 0;
329 : }
330 60 : myNet.getThreadPool().add(new RONet::BulkmodeTask(false), workerIndex);
331 : lastOrigin = c->origin;
332 60 : myNet.getThreadPool().add(new RoutingTask(*this, c, begin, linkFlow), workerIndex);
333 : } else {
334 0 : myNet.getThreadPool().add(new RoutingTask(*this, c, begin, linkFlow), workerIndex);
335 : }
336 60 : continue;
337 : }
338 : #else
339 : UNUSED_PARAMETER(workerIndex);
340 : #endif
341 40800 : if (lastOrigin != c->origin) {
342 11960 : myRouter.setBulkMode(false);
343 : lastOrigin = c->origin;
344 : }
345 40800 : computePath(c, begin, linkFlow, &myRouter, true);
346 : }
347 : #ifdef HAVE_FOX
348 4230 : if (myNet.getThreadPool().size() > 0) {
349 40 : myNet.getThreadPool().waitAll();
350 : }
351 : #endif
352 45090 : for (std::vector<ODCell*>::const_iterator i = firstCell; i != lastCell; i++) {
353 40860 : ODCell* const c = *i;
354 40860 : const double linkFlow = c->vehicleNumber / numIter;
355 40860 : const SUMOTime begin = myAdditiveTraffic ? myBegin : c->begin;
356 40860 : const SUMOTime end = myAdditiveTraffic ? myEnd : c->end;
357 40860 : const double intervalLengthInHours = STEPS2TIME(end - begin) / 3600.;
358 40860 : const ConstROEdgeVector& edges = c->pathsVector.back()->getEdgeVector();
359 283024 : for (ConstROEdgeVector::const_iterator e = edges.begin(); e != edges.end(); e++) {
360 242164 : ROMAEdge* edge = static_cast<ROMAEdge*>(myNet.getEdge((*e)->getID()));
361 242164 : const double newFlow = edge->getFlow(STEPS2TIME(begin)) + linkFlow;
362 242164 : edge->setFlow(STEPS2TIME(begin), STEPS2TIME(end), newFlow);
363 242164 : double travelTime = capacityConstraintFunction(edge, newFlow / intervalLengthInHours);
364 242164 : if (lastBegin >= 0 && myAdaptionFactor > 0.) {
365 : if (loadedTravelTimes.count(edge) != 0) {
366 0 : travelTime = loadedTravelTimes[edge] * myAdaptionFactor + (1. - myAdaptionFactor) * travelTime;
367 : } else {
368 55200 : travelTime = edge->getTravelTime(myDefaultVehicle, STEPS2TIME(lastBegin)) * myAdaptionFactor + (1. - myAdaptionFactor) * travelTime;
369 : }
370 : }
371 242164 : edge->addTravelTime(travelTime, STEPS2TIME(begin), STEPS2TIME(end));
372 : }
373 : }
374 4230 : myRouter.reset(myDefaultVehicle);
375 : }
376 210 : writeInterval(intervalStart, (*(lastCell - 1))->end);
377 : lastBegin = intervalStart;
378 : }
379 69 : }
380 :
381 :
382 : void
383 4 : ROMAAssignments::sue(const int maxOuterIteration, const int maxInnerIteration, const int kPaths, const double penalty, const double tolerance, const std::string /* routeChoiceMethod */) {
384 4 : getKPaths(kPaths, penalty);
385 : std::map<const SUMOTime, SUMOTime> intervals;
386 4 : if (myAdditiveTraffic) {
387 0 : intervals[myBegin] = myEnd;
388 : } else {
389 8 : for (const ODCell* const c : myMatrix.getCells()) {
390 4 : intervals[c->begin] = c->end;
391 : }
392 : }
393 5 : for (int outer = 0; outer < maxOuterIteration; outer++) {
394 11 : for (int inner = 0; inner < maxInnerIteration; inner++) {
395 22 : for (const ODCell* const c : myMatrix.getCells()) {
396 11 : const SUMOTime begin = myAdditiveTraffic ? myBegin : c->begin;
397 11 : const SUMOTime end = myAdditiveTraffic ? myEnd : c->end;
398 : // update path cost
399 29 : for (std::vector<RORoute*>::const_iterator j = c->pathsVector.begin(); j != c->pathsVector.end(); ++j) {
400 18 : RORoute* r = *j;
401 18 : r->setCosts(myRouter.recomputeCosts(r->getEdgeVector(), myDefaultVehicle, 0));
402 : // std::cout << std::setprecision(20) << r->getID() << ":" << r->getCosts() << std::endl;
403 : }
404 : // calculate route utilities and probabilities
405 22 : RouteCostCalculator<RORoute, ROEdge, ROVehicle>::getCalculator().calculateProbabilities(c->pathsVector, myDefaultVehicle, 0);
406 : // calculate route flows
407 29 : for (std::vector<RORoute*>::const_iterator j = c->pathsVector.begin(); j != c->pathsVector.end(); ++j) {
408 18 : RORoute* r = *j;
409 18 : const double pathFlow = r->getProbability() * c->vehicleNumber;
410 : // assign edge flow deltas
411 201 : for (ConstROEdgeVector::const_iterator e = r->getEdgeVector().begin(); e != r->getEdgeVector().end(); e++) {
412 183 : ROMAEdge* edge = static_cast<ROMAEdge*>(myNet.getEdge((*e)->getID()));
413 183 : edge->setHelpFlow(STEPS2TIME(begin), STEPS2TIME(end), edge->getHelpFlow(STEPS2TIME(begin)) + pathFlow);
414 : }
415 : }
416 : }
417 : // calculate new edge flows and check for stability
418 : int unstableEdges = 0;
419 22 : for (const auto& it : intervals) {
420 11 : const double intervalLengthInHours = STEPS2TIME(it.second - it.first) / 3600.;
421 11 : const double intBegin = STEPS2TIME(it.first);
422 11 : const double intEnd = STEPS2TIME(it.second);
423 447 : for (ROEdge* e : ROEdge::getAllEdges()) {
424 : ROMAEdge* const edge = static_cast<ROMAEdge*>(e);
425 : const double oldFlow = edge->getFlow(intBegin);
426 : double newFlow = oldFlow;
427 436 : if (inner == 0 && outer == 0) {
428 149 : newFlow += edge->getHelpFlow(intBegin);
429 : } else {
430 287 : newFlow += (edge->getHelpFlow(intBegin) - oldFlow) / (inner + 1);
431 : }
432 : // if not lohse:
433 436 : if (newFlow > 0.) {
434 137 : if (fabs(newFlow - oldFlow) / newFlow > tolerance) {
435 60 : unstableEdges++;
436 : }
437 299 : } else if (newFlow == 0.) {
438 299 : if (oldFlow != 0. && (fabs(newFlow - oldFlow) / oldFlow > tolerance)) {
439 0 : unstableEdges++;
440 : }
441 : } else { // newFlow < 0.
442 0 : unstableEdges++;
443 : newFlow = 0.;
444 : }
445 : edge->setFlow(intBegin, intEnd, newFlow);
446 436 : const double travelTime = capacityConstraintFunction(edge, newFlow / intervalLengthInHours);
447 436 : edge->addTravelTime(travelTime, intBegin, intEnd);
448 : edge->setHelpFlow(intBegin, intEnd, 0.);
449 : }
450 : }
451 : // if stable break
452 11 : if (unstableEdges == 0) {
453 : break;
454 : }
455 : // additional stability check from python script: if notstable < math.ceil(net.geteffEdgeCounts()*0.005) or notstable < 3: stable = True
456 : }
457 : // check for a new route, if none available, break
458 : // several modifications about when a route is new and when to break are in the original script
459 : bool newRoute = false;
460 10 : for (ODCell* const c : myMatrix.getCells()) {
461 5 : newRoute |= !computePath(c).empty();
462 : }
463 5 : if (!newRoute) {
464 : break;
465 : }
466 : }
467 : // final round of assignment
468 8 : for (const ODCell* const c : myMatrix.getCells()) {
469 : // update path cost
470 11 : for (std::vector<RORoute*>::const_iterator j = c->pathsVector.begin(); j != c->pathsVector.end(); ++j) {
471 7 : RORoute* r = *j;
472 7 : r->setCosts(myRouter.recomputeCosts(r->getEdgeVector(), myDefaultVehicle, 0));
473 : }
474 : // calculate route utilities and probabilities
475 8 : RouteCostCalculator<RORoute, ROEdge, ROVehicle>::getCalculator().calculateProbabilities(c->pathsVector, myDefaultVehicle, 0);
476 : // calculate route flows
477 11 : for (std::vector<RORoute*>::const_iterator j = c->pathsVector.begin(); j != c->pathsVector.end(); ++j) {
478 7 : RORoute* r = *j;
479 7 : r->setProbability(r->getProbability() * c->vehicleNumber);
480 : }
481 : }
482 8 : for (const auto& it : intervals) {
483 4 : writeInterval(it.first, it.second);
484 : }
485 4 : }
486 :
487 :
488 : double
489 0 : ROMAAssignments::getPenalizedEffort(const ROEdge* const e, const ROVehicle* const v, double t) {
490 : const std::map<const ROEdge* const, double>::const_iterator i = myPenalties.find(e);
491 0 : return i == myPenalties.end() ? e->getEffort(v, t) : e->getEffort(v, t) + i->second;
492 : }
493 :
494 :
495 : double
496 236 : ROMAAssignments::getPenalizedTT(const ROEdge* const e, const ROVehicle* const v, double t) {
497 : const std::map<const ROEdge* const, double>::const_iterator i = myPenalties.find(e);
498 236 : return i == myPenalties.end() ? e->getTravelTime(v, t) : e->getTravelTime(v, t) + i->second;
499 : }
500 :
501 :
502 : double
503 0 : ROMAAssignments::getTravelTime(const ROEdge* const e, const ROVehicle* const v, double t) {
504 0 : return e->getTravelTime(v, t);
505 : }
506 :
507 :
508 : #ifdef HAVE_FOX
509 : // ---------------------------------------------------------------------------
510 : // ROMAAssignments::RoutingTask-methods
511 : // ---------------------------------------------------------------------------
512 : void
513 60 : ROMAAssignments::RoutingTask::run(MFXWorkerThread* context) {
514 120 : myAssign.computePath(myCell, myBegin, myLinkFlow, &static_cast<RONet::WorkerThread*>(context)->getVehicleRouter(SVC_IGNORING), mySetBulkMode);
515 60 : }
516 : #endif
|
