Eclipse SUMO - Simulation of Urban MObility
Loading...
Searching...
No Matches
Circuit.cpp
Go to the documentation of this file.
1/****************************************************************************/
2// Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
3// Copyright (C) 2001-2023 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/****************************************************************************/
24// Representation of electric circuit of overhead wires
25/****************************************************************************/
26#include <config.h>
27
28#include <cfloat>
29#include <cstdlib>
30#include <iostream>
31#include <ctime>
32#include <mutex>
35#include <microsim/MSGlobals.h>
36#include "Element.h"
37#include "Node.h"
38#include "Circuit.h"
39
40static std::mutex circuit_lock;
41
42Node* Circuit::addNode(std::string name) {
43 if (getNode(name) != nullptr) {
44 WRITE_ERRORF(TL("The node: '%' already exists."), name);
45 return nullptr;
46 }
47
48 if (nodes->size() == 0) {
49 lastId = -1;
50 }
51 Node* tNode = new Node(name, this->lastId);
52 if (lastId == -1) {
53 tNode->setGround(true);
54 }
55 this->lastId++;
56 circuit_lock.lock();
57 this->nodes->push_back(tNode);
58 circuit_lock.unlock();
59 return tNode;
60}
61
63 circuit_lock.lock();
64 this->nodes->erase(std::remove(this->nodes->begin(), this->nodes->end(), node), this->nodes->end());
65 circuit_lock.unlock();
66}
67
68double Circuit::getCurrent(std::string name) {
69 Element* tElement = getElement(name);
70 if (tElement == nullptr) {
71 return DBL_MAX;
72 }
73 return tElement->getCurrent();
74}
75
76double Circuit::getVoltage(std::string name) {
77 Element* tElement = getElement(name);
78 if (tElement == nullptr) {
79 Node* node = getNode(name);
80 if (node != nullptr) {
81 return node->getVoltage();
82 } else {
83 return DBL_MAX;
84 }
85 } else {
86 return tElement->getVoltage();
87 }
88}
89
90double Circuit::getResistance(std::string name) {
91 Element* tElement = getElement(name);
92 if (tElement == nullptr) {
93 return -1;
94 }
95 return tElement->getResistance();
96}
97
98Node* Circuit::getNode(std::string name) {
99 for (Node* const node : *nodes) {
100 if (node->getName() == name) {
101 return node;
102 }
103 }
104 return nullptr;
105}
106
108 for (Node* const node : *nodes) {
109 if (node->getId() == id) {
110 return node;
111 }
112 }
113 return nullptr;
114}
115
116Element* Circuit::getElement(std::string name) {
117 for (Element* const el : *elements) {
118 if (el->getName() == name) {
119 return el;
120 }
121 }
122 for (Element* const voltageSource : *voltageSources) {
123 if (voltageSource->getName() == name) {
124 return voltageSource;
125 }
126 }
127 return nullptr;
128}
129
131 for (Element* const el : *elements) {
132 if (el->getId() == id) {
133 return el;
134 }
135 }
136 return getVoltageSource(id);
137}
138
140 for (Element* const voltageSource : *voltageSources) {
141 if (voltageSource->getId() == id) {
142 return voltageSource;
143 }
144 }
145 return nullptr;
146}
147
149 double power = 0;
150 for (Element* const voltageSource : *voltageSources) {
151 power += voltageSource->getPower();
152 }
153 return power;
154}
155
157 double current = 0;
158 for (Element* const voltageSource : *voltageSources) {
159 current += voltageSource->getCurrent();
160 }
161 return current;
162}
163
164// RICE_CHECK: Locking removed?
165std::string& Circuit::getCurrentsOfCircuitSource(std::string& currents) {
166 //circuit_lock.lock();
167 currents.clear();
168 for (Element* const voltageSource : *voltageSources) {
169 currents += toString(voltageSource->getCurrent(), 4) + " ";
170 }
171 if (!currents.empty()) {
172 currents.pop_back();
173 }
174 //circuit_lock.unlock();
175 return currents;
176}
177
178std::vector<Element*>* Circuit::getCurrentSources() {
179 std::vector<Element*>* vsources = new std::vector<Element*>(0);
180 for (Element* const el : *elements) {
182 //if ((*it)->getType() == Element::ElementType::CURRENT_SOURCE_traction_wire && !isnan((*it)->getPowerWanted())) {
183 vsources->push_back(el);
184 }
185 }
186 return vsources;
187}
188
190 circuit_lock.lock();
191}
192
194 circuit_lock.unlock();
195}
196
197#ifdef HAVE_EIGEN
198void Circuit::removeColumn(Eigen::MatrixXd& matrix, int colToRemove) {
199 const int numRows = (int)matrix.rows();
200 const int numCols = (int)matrix.cols() - 1;
201
202 if (colToRemove < numCols) {
203 matrix.block(0, colToRemove, numRows, numCols - colToRemove) = matrix.rightCols(numCols - colToRemove);
204 }
205
206 matrix.conservativeResize(numRows, numCols);
207}
208
209bool Circuit::solveEquationsNRmethod(double* eqn, double* vals, std::vector<int>* removable_ids) {
210 // removable_ids includes nodes with voltage source already
211 int numofcolumn = (int)voltageSources->size() + (int)nodes->size() - 1;
212 int numofeqs = numofcolumn - (int)removable_ids->size();
213
214 // map equations into matrix A
215 Eigen::MatrixXd A = Eigen::Map < Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> >(eqn, numofeqs, numofcolumn);
216
217 int id;
218 // remove removable columns of matrix A, i.e. remove equations corresponding to nodes with two resistors connected in series
219 // RICE_TODO auto for ?
220 for (std::vector<int>::reverse_iterator it = removable_ids->rbegin(); it != removable_ids->rend(); ++it) {
221 id = (*it >= 0 ? *it : -(*it));
222 removeColumn(A, id);
223 }
224
225 // detect number of column for each node
226 // in other words: detect elements of x to certain node
227 // in other words: assign number of column to the proper non removable node
228 int j = 0;
229 Element* tElem = nullptr;
230 Node* tNode = nullptr;
231 for (int i = 0; i < numofcolumn; i++) {
232 tNode = getNode(i);
233 if (tNode != nullptr) {
234 if (tNode->isRemovable()) {
235 tNode->setNumMatrixCol(-1);
236 continue;
237 } else {
238 if (j > numofeqs) {
239 WRITE_ERROR(TL("Index of renumbered node exceeded the reduced number of equations."));
240 break;
241 }
242 tNode->setNumMatrixCol(j);
243 j++;
244 continue;
245 }
246 } else {
247 tElem = getElement(i);
248 if (tElem != nullptr) {
249 if (j > numofeqs) {
250 WRITE_ERROR(TL("Index of renumbered element exceeded the reduced number of equations."));
251 break;
252 }
253 continue;
254 }
255 }
256 // tNode == nullptr && tElem == nullptr
257 WRITE_ERROR(TL("Structural error in reduced circuit matrix."));
258 }
259
260 // map 'vals' into vector b and initialize solution x
261 Eigen::Map<Eigen::VectorXd> b(vals, numofeqs);
262 Eigen::VectorXd x = A.colPivHouseholderQr().solve(b);
263
264 // initialize Jacobian matrix J and vector dx
265 Eigen::MatrixXd J = A;
266 Eigen::VectorXd dx;
267 // initialize progressively increasing maximal number of Newton-Rhapson iterations
268 int max_iter_of_NR = 10;
269 // number of tested values of alpha
270 int attemps = 0;
271 // value of scaling parameter alpha
272 double alpha = 1;
273 // the best (maximum) value of alpha that guarantees the existence of solution
274 alphaBest = 0;
275 // reason why is alpha not 1
277 // vector of alphas for that no solution has been found
278 std::vector<double> alpha_notSolution;
279 // initialize progressively decreasing tolerance for alpha
280 double alpha_res = 1e-2;
281
282 double currentSumActual = 0.0;
283 // solution x corresponding to the alphaBest
284 Eigen::VectorXd x_best = x;
285 bool x_best_exist = true;
286
287 if (x.maxCoeff() > 10e6 || x.minCoeff() < -10e6) {
288 WRITE_ERROR(TL("Initial solution x used during solving DC circuit is out of bounds.\n"));
289 }
290
291 // Search for the suitable scaling value alpha
292 while (true) {
293
294 ++attemps;
295 int iterNR = 0;
296 // run Newton-Raphson methods
297 while (true) {
298
299 // update right-hand side vector vals and Jacobian matrix J
300 // node's right-hand side set to zero
301 for (int i = 0; i < numofeqs - (int) voltageSources->size(); i++) {
302 vals[i] = 0;
303 }
304 J = A;
305
306 int i = 0;
307 for (auto& node : *nodes) {
308 if (node->isGround() || node->isRemovable() || node->getNumMatrixRow() == -2) {
309 continue;
310 }
311 if (node->getNumMatrixRow() != i) {
312 WRITE_ERROR(TL("wrongly assigned row of matrix A during solving the circuit"));
313 }
314 // TODO: Range-based loop
315 // loop over all node's elements
316 for (auto it_element = node->getElements()->begin(); it_element != node->getElements()->end(); it_element++) {
317 if ((*it_element)->getType() == Element::ElementType::CURRENT_SOURCE_traction_wire) {
318 // if the element is current source
319 if ((*it_element)->isEnabled()) {
320 double diff_voltage;
321 int PosNode_NumACol = (*it_element)->getPosNode()->getNumMatrixCol();
322 int NegNode_NumACol = (*it_element)->getNegNode()->getNumMatrixCol();
323 // compute voltage on current source
324 if (PosNode_NumACol == -1) {
325 // if the positive node is the ground => U = 0 - phi(NegNode)
326 diff_voltage = -x[NegNode_NumACol];
327 } else if (NegNode_NumACol == -1) {
328 // if the negative node is the ground => U = phi(PosNode) - 0
329 diff_voltage = x[PosNode_NumACol];
330 } else {
331 // U = phi(PosNode) - phi(NegNode)
332 diff_voltage = (x[PosNode_NumACol] - x[NegNode_NumACol]);
333 }
334
335 if ((*it_element)->getPosNode() == node) {
336 // the positive current (the element is consuming energy if powerWanted > 0) is flowing from the positive node (sign minus)
337 vals[i] -= alpha * (*it_element)->getPowerWanted() / diff_voltage;
338 (*it_element)->setCurrent(-alpha * (*it_element)->getPowerWanted() / diff_voltage);
339 if (PosNode_NumACol != -1) {
340 // -1* d_b/d_phiPos = -1* d(-alpha*P/(phiPos-phiNeg) )/d_phiPos = -1* (--alpha*P/(phiPos-phiNeg)^2 )
341 J(i, PosNode_NumACol) -= alpha * (*it_element)->getPowerWanted() / diff_voltage / diff_voltage;
342 }
343 if (NegNode_NumACol != -1) {
344 // -1* d_b/d_phiNeg = -1* d(-alpha*P/(phiPos-phiNeg) )/d_phiNeg = -1* (---alpha*P/(phiPos-phiNeg)^2 )
345 J(i, NegNode_NumACol) += alpha * (*it_element)->getPowerWanted() / diff_voltage / diff_voltage;
346 }
347 } else {
348 // the positive current (the element is consuming energy if powerWanted > 0) is flowing to the negative node (sign plus)
349 vals[i] += alpha * (*it_element)->getPowerWanted() / diff_voltage;
350 //Question: sign before alpha - or + during setting current?
351 //Answer: sign before alpha is minus since we assume positive powerWanted if the current element behaves as load
352 // (*it_element)->setCurrent(-alpha * (*it_element)->getPowerWanted() / diff_voltage);
353 // Note: we should never reach this part of code since the authors assumes the negataive node of current source as the ground node
354 WRITE_WARNING(TL("The negative node of current source is not the groud."))
355 if (PosNode_NumACol != -1) {
356 // -1* d_b/d_phiPos = -1* d(alpha*P/(phiPos-phiNeg) )/d_phiPos = -1* (-alpha*P/(phiPos-phiNeg)^2 )
357 J(i, PosNode_NumACol) += alpha * (*it_element)->getPowerWanted() / diff_voltage / diff_voltage;
358 }
359 if (NegNode_NumACol != -1) {
360 // -1* d_b/d_phiNeg = -1* d(alpha*P/(phiPos-phiNeg) )/d_phiNeg = -1* (--alpha*P/(phiPos-phiNeg)^2 )
361 J(i, NegNode_NumACol) -= alpha * (*it_element)->getPowerWanted() / diff_voltage / diff_voltage;
362 }
363 }
364 }
365 }
366 }
367 i++;
368 }
369
370
371 // RICE_CHECK @20210409 This had to be merged into the master/main manually.
372 // Sum of currents going through the all voltage sources
373 // the sum is over all nodes, but the nonzero nodes are only those neigboring with current sources,
374 // so the sum is negative sum of currents through/from current sources representing trolleybusess
375 currentSumActual = 0;
376 for (i = 0; i < numofeqs - (int)voltageSources->size(); i++) {
377 currentSumActual -= vals[i];
378 }
379 // RICE_TODO @20210409 This epsilon should be specified somewhere as a constant. Or should be a parameter.
380 if ((A * x - b).norm() < 1e-6) {
381 //current limits
382 if (currentSumActual > getCurrentLimit() && MSGlobals::gOverheadWireCurrentLimits) {
384 alpha_notSolution.push_back(alpha);
385 if (x_best_exist) {
386 x = x_best;
387 }
388 break;
389 }
390 //voltage limits 70% - 120% of nominal voltage
391 // RICE_TODO @20210409 Again, these limits should be parametrised.
392 if (x.maxCoeff() > voltageSources->front()->getVoltage() * 1.2 || x.minCoeff() < voltageSources->front()->getVoltage() * 0.7) {
394 alpha_notSolution.push_back(alpha);
395 if (x_best_exist) {
396 x = x_best;
397 }
398 break;
399 }
400
401 alphaBest = alpha;
402 x_best = x;
403 x_best_exist = true;
404 break;
405 } else if (iterNR == max_iter_of_NR) {
407 alpha_notSolution.push_back(alpha);
408 if (x_best_exist) {
409 x = x_best;
410 }
411 break;
412 }
413
414 // Newton=Rhapson iteration
415 dx = -J.colPivHouseholderQr().solve(A * x - b);
416 x = x + dx;
417 ++iterNR;
418 }
419
420 if (alpha_notSolution.empty()) {
421 // no alpha without solution is in the alpha_notSolution, so the solving procedure is terminating
422 break;
423 }
424
425 if ((alpha_notSolution.back() - alphaBest) < alpha_res) {
426 max_iter_of_NR = 2 * max_iter_of_NR;
427 // RICE_TODO @20210409 Why division by 10?
428 // it follows Sevcik, Jakub, et al. "Solvability of the Power Flow Problem in DC Overhead Wire Circuit Modeling." Applications of Mathematics (2021): 1-19.
429 // see Alg 2 (progressive decrease of optimality tolerance)
430 alpha_res = alpha_res / 10;
431 // RICE_TODO @20210409 This epsilon should be specified somewhere as a constant. Or should be a parameter.
432 if (alpha_res < 5e-5) {
433 break;
434 }
435 alpha = alpha_notSolution.back();
436 alpha_notSolution.pop_back();
437 continue;
438 }
439
440 alpha = alphaBest + 0.5 * (alpha_notSolution.back() - alphaBest);
441 }
442
443 // vals is pointer to memory and we use it now for saving solution x_best instead of right-hand side b
444 for (int i = 0; i < numofeqs; i++) {
445 vals[i] = x_best[i];
446 }
447
448 // RICE_TODO: Describe what is hapenning here.
449 // we take x_best and alphaBest and update current values in current sources in order to be in agreement with the solution
450 int i = 0;
451 for (auto& node : *nodes) {
452 if (node->isGround() || node->isRemovable() || node->getNumMatrixRow() == -2) {
453 continue;
454 }
455 if (node->getNumMatrixRow() != i) {
456 WRITE_ERROR(TL("wrongly assigned row of matrix A during solving the circuit"));
457 }
458 for (auto it_element = node->getElements()->begin(); it_element != node->getElements()->end(); it_element++) {
459 if ((*it_element)->getType() == Element::ElementType::CURRENT_SOURCE_traction_wire) {
460 if ((*it_element)->isEnabled()) {
461 double diff_voltage;
462 int PosNode_NumACol = (*it_element)->getPosNode()->getNumMatrixCol();
463 int NegNode_NumACol = (*it_element)->getNegNode()->getNumMatrixCol();
464 if (PosNode_NumACol == -1) {
465 diff_voltage = -x_best[NegNode_NumACol];
466 } else if (NegNode_NumACol == -1) {
467 diff_voltage = x_best[PosNode_NumACol];
468 } else {
469 diff_voltage = (x_best[PosNode_NumACol] - x_best[NegNode_NumACol]);
470 }
471
472 if ((*it_element)->getPosNode() == node) {
473 (*it_element)->setCurrent(-alphaBest * (*it_element)->getPowerWanted() / diff_voltage);
474 } else {
475 //Question: sign before alpha - or + during setting current?
476 //Answer: sign before alpha is minus since we assume positive powerWanted if the current element behaves as load
477 // (*it_element)->setCurrent(-alphaBest * (*it_element)->getPowerWanted() / diff_voltage);
478 // Note: we should never reach this part of code since the authors assumes the negataive node of current source as the ground node
479 WRITE_WARNING(TL("The negative node of current source is not the groud."))
480 }
481 }
482 }
483 }
484 i++;
485 }
486
487 return true;
488}
489#endif
490
491void Circuit::deployResults(double* vals, std::vector<int>* removable_ids) {
492 // vals are the solution x
493
494 int numofcolumn = (int)voltageSources->size() + (int)nodes->size() - 1;
495 int numofeqs = numofcolumn - (int)removable_ids->size();
496
497 //loop over non-removable nodes: we assign the computed voltage to the non-removables nodes
498 int j = 0;
499 Element* tElem = nullptr;
500 Node* tNode = nullptr;
501 for (int i = 0; i < numofcolumn; i++) {
502 tNode = getNode(i);
503 if (tNode != nullptr)
504 if (tNode->isRemovable()) {
505 continue;
506 } else {
507 if (j > numofeqs) {
508 WRITE_ERROR(TL("Results deployment during circuit evaluation was unsuccessful."));
509 break;
510 }
511 tNode->setVoltage(vals[j]);
512 j++;
513 continue;
514 } else {
515 tElem = getElement(i);
516 if (tElem != nullptr) {
517 if (j > numofeqs) {
518 WRITE_ERROR(TL("Results deployment during circuit evaluation was unsuccessful."));
519 break;
520 }
521 // tElem should be voltage source - the current through voltage source is computed in a loop below
522 // if tElem is current source (JS thinks that no current source's id <= numofeqs), the current is already assign at the end of solveEquationsNRmethod method
523 continue;
524 }
525 }
526 WRITE_ERROR(TL("Results deployment during circuit evaluation was unsuccessful."));
527 }
528
529 Element* el1 = nullptr;
530 Element* el2 = nullptr;
531 Node* nextNONremovableNode1 = nullptr;
532 Node* nextNONremovableNode2 = nullptr;
533 // interpolate result of voltage to removable nodes
534 for (Node* const node : *nodes) {
535 if (!node->isRemovable()) {
536 continue;
537 }
538 if (node->getElements()->size() != 2) {
539 continue;
540 }
541
542 el1 = node->getElements()->front();
543 el2 = node->getElements()->back();
544 nextNONremovableNode1 = el1->getTheOtherNode(node);
545 nextNONremovableNode2 = el2->getTheOtherNode(node);
546 double x = el1->getResistance();
547 double y = el2->getResistance();
548
549 while (nextNONremovableNode1->isRemovable()) {
550 el1 = nextNONremovableNode1->getAnOtherElement(el1);
551 x += el1->getResistance();
552 nextNONremovableNode1 = el1->getTheOtherNode(nextNONremovableNode1);
553 }
554
555 while (nextNONremovableNode2->isRemovable()) {
556 el2 = nextNONremovableNode2->getAnOtherElement(el2);
557 y += el2->getResistance();
558 nextNONremovableNode2 = el2->getTheOtherNode(nextNONremovableNode2);
559 }
560
561 x = x / (x + y);
562 y = ((1 - x) * nextNONremovableNode1->getVoltage()) + (x * nextNONremovableNode2->getVoltage());
563 node->setVoltage(((1 - x)*nextNONremovableNode1->getVoltage()) + (x * nextNONremovableNode2->getVoltage()));
564 node->setRemovability(false);
565 }
566
567 // Update the electric currents for voltage sources (based on Kirchhof's law: current out = current in)
568 for (Element* const voltageSource : *voltageSources) {
569 double currentSum = 0;
570 for (Element* const el : *voltageSource->getPosNode()->getElements()) {
571 // loop over all elements on PosNode excluding the actual voltage source it
572 if (el != voltageSource) {
573 //currentSum += el->getCurrent();
574 currentSum += (voltageSource->getPosNode()->getVoltage() - el->getTheOtherNode(voltageSource->getPosNode())->getVoltage()) / el->getResistance();
576 WRITE_WARNING(TL("Cannot assign unambigous electric current value to two voltage sources connected in parallel at the same node."));
577 }
578 }
579 }
580 voltageSource->setCurrent(currentSum);
581 }
582}
583
585 nodes = new std::vector<Node*>(0);
586 elements = new std::vector<Element*>(0);
587 voltageSources = new std::vector<Element*>(0);
588 lastId = 0;
589 iscleaned = true;
590 circuitCurrentLimit = INFINITY;
591}
592
593Circuit::Circuit(double currentLimit) {
594 nodes = new std::vector<Node*>(0);
595 elements = new std::vector<Element*>(0);
596 voltageSources = new std::vector<Element*>(0);
597 lastId = 0;
598 iscleaned = true;
599 circuitCurrentLimit = currentLimit;
600}
601
602#ifdef HAVE_EIGEN
603bool Circuit::_solveNRmethod() {
604 double* eqn = nullptr;
605 double* vals = nullptr;
606 std::vector<int> removable_ids;
607
608 detectRemovableNodes(&removable_ids);
609 createEquationsNRmethod(eqn, vals, &removable_ids);
610 if (!solveEquationsNRmethod(eqn, vals, &removable_ids)) {
611 return false;
612 }
613 // vals are now the solution x of the circuit
614 deployResults(vals, &removable_ids);
615
616 delete eqn;
617 delete vals;
618 return true;
619}
620
621bool Circuit::solve() {
622 if (!iscleaned) {
623 cleanUpSP();
624 }
625 return this->_solveNRmethod();
626}
627
628bool Circuit::createEquationsNRmethod(double*& eqs, double*& vals, std::vector<int>* removable_ids) {
629 // removable_ids does not include nodes with voltage source yet
630
631 // number of voltage sources + nodes without the ground node
632 int n = (int)(voltageSources->size() + nodes->size() - 1);
633 // number of equations
634 // assumption: each voltage source has different positive node and common ground node,
635 // i.e. any node excluding the ground node is connected to 0 or 1 voltage source
636 int m = n - (int)(removable_ids->size() + voltageSources->size());
637
638 // allocate and initialize zero matrix eqs and vector vals
639 eqs = new double[m * n];
640 vals = new double[m];
641
642 for (int i = 0; i < m; i++) {
643 vals[i] = 0;
644 for (int j = 0; j < n; j++) {
645 eqs[i * n + j] = 0;
646 }
647 }
648
649 // loop over all nodes
650 int i = 0;
651 for (std::vector<Node*>::iterator it = nodes->begin(); it != nodes->end(); it++) {
652 if ((*it)->isGround() || (*it)->isRemovable()) {
653 // if the node is grounded or is removable set the corresponding number of row in matrix to -1 (no equation in eqs)
654 (*it)->setNumMatrixRow(-1);
655 continue;
656 }
657 assert(i < m);
658 // constitute the equation corresponding to node it, add all passed voltage source elements into removable_ids
659 bool noVoltageSource = createEquationNRmethod((*it), (eqs + n * i), vals[i], removable_ids);
660 // if the node it has element of type "voltage source" we do not use the equation, because some value of current throw the voltage source can be always find
661 if (noVoltageSource) {
662 (*it)->setNumMatrixRow(i);
663 i++;
664 } else {
665 (*it)->setNumMatrixRow(-2);
666 vals[i] = 0;
667 for (int j = 0; j < n; j++) {
668 eqs[n * i + j] = 0;
669 }
670 }
671 }
672
673 // removable_ids includes nodes with voltage source already
674 std::sort(removable_ids->begin(), removable_ids->end(), std::less<int>());
675
676
677 for (std::vector<Element*>::iterator it = voltageSources->begin(); it != voltageSources->end(); it++) {
678 assert(i < m);
679 createEquation((*it), (eqs + n * i), vals[i]);
680 i++;
681 }
682
683 return true;
684}
685
686bool Circuit::createEquation(Element* vsource, double* eqn, double& val) {
687 if (!vsource->getPosNode()->isGround()) {
688 eqn[vsource->getPosNode()->getId()] = 1;
689 }
690 if (!vsource->getNegNode()->isGround()) {
691 eqn[vsource->getNegNode()->getId()] = -1;
692 }
693 if (vsource->isEnabled()) {
694 val = vsource->getVoltage();
695 } else {
696 val = 0;
697 }
698 return true;
699}
700
701bool Circuit::createEquationNRmethod(Node* node, double* eqn, double& val, std::vector<int>* removable_ids) {
702 // loop over all elements connected to the node
703 for (std::vector<Element*>::iterator it = node->getElements()->begin(); it != node->getElements()->end(); it++) {
704 double x;
705 switch ((*it)->getType()) {
707 if ((*it)->isEnabled()) {
708 x = (*it)->getResistance();
709 // go through all neigboring removable nodes and sum resistance of resistors in the serial branch
710 Node* nextNONremovableNode = (*it)->getTheOtherNode(node);
711 Element* nextSerialResistor = *it;
712 while (nextNONremovableNode->isRemovable()) {
713 nextSerialResistor = nextNONremovableNode->getAnOtherElement(nextSerialResistor);
714 x += nextSerialResistor->getResistance();
715 nextNONremovableNode = nextSerialResistor->getTheOtherNode(nextNONremovableNode);
716 }
717 // compute inverse value and place/add this value at proper places in eqn
718 x = 1 / x;
719 eqn[node->getId()] += x;
720
721 if (!nextNONremovableNode->isGround()) {
722 eqn[nextNONremovableNode->getId()] -= x;
723 }
724 }
725 break;
727 if ((*it)->isEnabled()) {
728 // initialize current in current source
729 if ((*it)->getPosNode() == node) {
730 x = -(*it)->getPowerWanted() / voltageSources->front()->getVoltage();
731 } else {
732 x = (*it)->getPowerWanted() / voltageSources->front()->getVoltage();
733 }
734 } else {
735 x = 0;
736 }
737 val += x;
738 break;
740 if ((*it)->getPosNode() == node) {
741 x = -1;
742 } else {
743 x = 1;
744 }
745 eqn[(*it)->getId()] += x;
746 // equations with voltage source can be ignored, because some value of current throw the voltage source can be always find
747 removable_ids->push_back((*it)->getId());
748 return false;
749 break;
751 return false;
752 break;
753 }
754 }
755 return true;
756}
757#endif
758
764void Circuit::detectRemovableNodes(std::vector<int>* removable_ids) {
765 // loop over all nodes in the circuit
766 for (std::vector<Node*>::iterator it = nodes->begin(); it != nodes->end(); it++) {
767 // if the node is connected to two elements and is not the ground
768 if ((*it)->getElements()->size() == 2 && !(*it)->isGround()) {
769 // set such node defaultly as removable. But check if the two elements are both resistors
770 (*it)->setRemovability(true);
771 for (std::vector<Element*>::iterator it2 = (*it)->getElements()->begin(); it2 != (*it)->getElements()->end(); it2++) {
772 if ((*it2)->getType() != Element::ElementType::RESISTOR_traction_wire) {
773 (*it)->setRemovability(false);
774 break;
775 }
776 }
777 if ((*it)->isRemovable()) {
778 //if the node is removeable add pointer into the vector of removeblas nodes
779 removable_ids->push_back((*it)->getId());
780 }
781 } else {
782 (*it)->setRemovability(false);
783 }
784 }
785 // sort the vector of removable ids
786 std::sort(removable_ids->begin(), removable_ids->end(), std::less<int>());
787 return;
788}
789
790Element* Circuit::addElement(std::string name, double value, Node* pNode, Node* nNode, Element::ElementType et) {
791 // RICE_CHECK: This seems to be a bit of work in progress, is it final?
792 // if ((et == Element::ElementType::RESISTOR_traction_wire && value <= 0) || et == Element::ElementType::ERROR_traction_wire) {
793 if (et == Element::ElementType::RESISTOR_traction_wire && value <= 1e-6) {
794 //due to numeric problems
795 // RICE_TODO @20210409 This epsilon should be specified somewhere as a constant. Or should be a parameter.
796 if (value > -1e-6) {
797 value = 1e-6;
798 WRITE_WARNING(TL("Trying to add resistor element into the overhead wire circuit with resistance < 1e-6. "))
799 } else {
800 WRITE_ERROR(TL("Trying to add resistor element into the overhead wire circuit with resistance < 0. "))
801 return nullptr;
802 }
803 }
804
805 Element* e = getElement(name);
806
807 if (e != nullptr) {
808 //WRITE_ERRORF(TL("The element: '%' already exists."), name);
809 std::cout << "The element: '" + name + "' already exists.";
810 return nullptr;
811 }
812
813 e = new Element(name, et, value);
815 e->setId(lastId);
816 lastId++;
817 circuit_lock.lock();
818 this->voltageSources->push_back(e);
819 circuit_lock.unlock();
820 } else {
821 circuit_lock.lock();
822 this->elements->push_back(e);
823 circuit_lock.unlock();
824 }
825
826 e->setPosNode(pNode);
827 e->setNegNode(nNode);
828
829 pNode->addElement(e);
830 nNode->addElement(e);
831 return e;
832}
833
835 element->getPosNode()->eraseElement(element);
836 element->getNegNode()->eraseElement(element);
837 circuit_lock.lock();
838 this->elements->erase(std::remove(this->elements->begin(), this->elements->end(), element), this->elements->end());
839 circuit_lock.unlock();
840}
841
842void Circuit::replaceAndDeleteNode(Node* unusedNode, Node* newNode) {
843 //replace element node if it is unusedNode
844 for (auto& voltageSource : *voltageSources) {
845 if (voltageSource->getNegNode() == unusedNode) {
846 voltageSource->setNegNode(newNode);
847 newNode->eraseElement(voltageSource);
848 newNode->addElement(voltageSource);
849 }
850 if (voltageSource->getPosNode() == unusedNode) {
851 voltageSource->setPosNode(newNode);
852 newNode->eraseElement(voltageSource);
853 newNode->addElement(voltageSource);
854 }
855 }
856 for (auto& element : *elements) {
857 if (element->getNegNode() == unusedNode) {
858 element->setNegNode(newNode);
859 newNode->eraseElement(element);
860 newNode->addElement(element);
861 }
862 if (element->getPosNode() == unusedNode) {
863 element->setPosNode(newNode);
864 newNode->eraseElement(element);
865 newNode->addElement(element);
866 }
867 }
868
869 //erase unusedNode from nodes vector
870 this->eraseNode(unusedNode);
871
872 //modify id of other elements and nodes
873 int modLastId = this->getLastId() - 1;
874 if (unusedNode->getId() != modLastId) {
875 Node* node_last = this->getNode(modLastId);
876 if (node_last != nullptr) {
877 node_last->setId(unusedNode->getId());
878 } else {
879 Element* elem_last = this->getVoltageSource(modLastId);
880 if (elem_last != nullptr) {
881 elem_last->setId(unusedNode->getId());
882 } else {
883 WRITE_ERROR(TL("The element or node with the last Id was not found in the circuit!"));
884 }
885 }
886 }
887
888 this->descreaseLastId();
889 delete unusedNode;
890}
891
893 for (std::vector<Element*>::iterator it = elements->begin(); it != elements->end(); it++) {
894 if ((*it)->getType() != Element::ElementType::RESISTOR_traction_wire) {
895 (*it)->setEnabled(true);
896 }
897 }
898
899 for (std::vector<Element*>::iterator it = voltageSources->begin(); it != voltageSources->end(); it++) {
900 (*it)->setEnabled(true);
901 }
902 this->iscleaned = true;
903}
904
905bool Circuit::checkCircuit(std::string substationId) {
906 // check empty nodes
907 for (std::vector<Node*>::iterator it = nodes->begin(); it != nodes->end(); it++) {
908 if ((*it)->getNumOfElements() < 2) {
909 //cout << "WARNING: Node [" << (*it)->getName() << "] is connected to less than two elements, please enter other elements.\n";
910 if ((*it)->getNumOfElements() < 1) {
911 return false;
912 }
913 }
914 }
915 // check voltage sources
916 for (std::vector<Element*>::iterator it = voltageSources->begin(); it != voltageSources->end(); it++) {
917 if ((*it)->getPosNode() == nullptr || (*it)->getNegNode() == nullptr) {
918 //cout << "ERROR: Voltage Source [" << (*it)->getName() << "] is connected to less than two nodes, please enter the other end.\n";
919 WRITE_ERRORF(TL("Circuit Voltage Source '%' is connected to less than two nodes, please adjust the definition of the section (with substation '%')."), (*it)->getName(), substationId);
920 return false;
921 }
922 }
923 // check other elements
924 for (std::vector<Element*>::iterator it = elements->begin(); it != elements->end(); it++) {
925 if ((*it)->getPosNode() == nullptr || (*it)->getNegNode() == nullptr) {
926 //cout << "ERROR: Element [" << (*it)->getName() << "] is connected to less than two nodes, please enter the other end.\n";
927 WRITE_ERRORF(TL("Circuit Element '%' is connected to less than two nodes, please adjust the definition of the section (with substation '%')."), (*it)->getName(), substationId);
928 return false;
929 }
930 }
931
932 // check connectivity
933 int num = (int)nodes->size() + getNumVoltageSources() - 1;
934 bool* nodesVisited = new bool[num];
935 for (int i = 0; i < num; i++) {
936 nodesVisited[i] = false;
937 }
938 // TODO: Probably unused
939 // int id = -1;
940 if (!getNode(-1)->isGround()) {
941 //cout << "ERROR: Node id -1 is not the ground \n";
942 WRITE_ERRORF(TL("Circuit Node with id '-1' is not the grounded, please adjust the definition of the section (with substation '%')."), substationId);
943 }
944 std::vector<Node*>* queue = new std::vector<Node*>(0);
945 Node* node = nullptr;
946 Node* neigboringNode = nullptr;
947 //start with (voltageSources->front()->getPosNode())
948 nodesVisited[voltageSources->front()->getId()] = 1;
949 node = voltageSources->front()->getPosNode();
950 queue->push_back(node);
951
952 while (!queue->empty()) {
953 node = queue->back();
954 queue->pop_back();
955 if (!nodesVisited[node->getId()]) {
956 nodesVisited[node->getId()] = true;
957 for (auto it = node->getElements()->begin(); it != node->getElements()->end(); it++) {
958 neigboringNode = (*it)->getTheOtherNode(node);
959 if (!neigboringNode->isGround()) {
960 queue->push_back(neigboringNode);
961 } else if ((*it)->getType() == Element::ElementType::VOLTAGE_SOURCE_traction_wire) {
963 nodesVisited[(*it)->getId()] = 1;
964 } else if ((*it)->getType() == Element::ElementType::RESISTOR_traction_wire) {
965 //cout << "ERROR: The resistor type connects the ground \n";
966 WRITE_ERRORF(TL("A Circuit Resistor Element connects the ground, please adjust the definition of the section (with substation '%')."), substationId);
967 }
968 }
969 }
970 }
971
972 for (int i = 0; i < num; i++) {
973 if (nodesVisited[i] == 0) {
974 //cout << "ERROR: Node or voltage source with id " << (i) << " has been not visited during checking of the circuit => Disconnectivity of the circuit. \n";
975 WRITE_WARNINGF(TL("Circuit Node or Voltage Source with internal id '%' has been not visited during checking of the circuit. The circuit is disconnected, please adjust the definition of the section (with substation '%')."), toString(i), substationId);
976 }
977 }
978
979 return true;
980}
981
983 return (int) voltageSources->size();
984}
static std::mutex circuit_lock
Definition Circuit.cpp:40
#define WRITE_WARNINGF(...)
Definition MsgHandler.h:271
#define WRITE_ERRORF(...)
Definition MsgHandler.h:280
#define WRITE_ERROR(msg)
Definition MsgHandler.h:279
#define WRITE_WARNING(msg)
Definition MsgHandler.h:270
#define TL(string)
Definition MsgHandler.h:287
std::string toString(const T &t, std::streamsize accuracy=gPrecision)
Definition ToString.h:46
std::vector< Node * > * nodes
Definition Circuit.h:66
std::vector< Element * > * getCurrentSources()
Definition Circuit.cpp:178
Node * addNode(std::string name)
Definition Circuit.cpp:42
double getVoltage(std::string name)
Definition Circuit.cpp:76
int getNumVoltageSources()
Definition Circuit.cpp:982
Element * addElement(std::string name, double value, Node *pNode, Node *nNode, Element::ElementType et)
Definition Circuit.cpp:790
int lastId
Definition Circuit.h:70
void eraseNode(Node *node)
Definition Circuit.cpp:62
void cleanUpSP()
Definition Circuit.cpp:892
double getCurrent(std::string name)
Definition Circuit.cpp:68
void unlock()
Definition Circuit.cpp:193
int getLastId()
Definition Circuit.h:240
Element * getElement(std::string name)
Definition Circuit.cpp:116
double circuitCurrentLimit
The electric current limit of the voltage sources.
Definition Circuit.h:74
double getTotalCurrentOfCircuitSources()
The sum of voltage source currents in the circuit.
Definition Circuit.cpp:156
void lock()
Definition Circuit.cpp:189
void detectRemovableNodes(std::vector< int > *removable_ids)
Definition Circuit.cpp:764
double getCurrentLimit()
@ brief Get the electric current limit of this circuit.
Definition Circuit.h:256
std::vector< Element * > * elements
Definition Circuit.h:67
void replaceAndDeleteNode(Node *unusedNode, Node *newNode)
Definition Circuit.cpp:842
Element * getVoltageSource(int id)
Definition Circuit.cpp:139
alphaFlag alphaReason
Definition Circuit.h:107
double getTotalPowerOfCircuitSources()
The sum of voltage source powers in the circuit.
Definition Circuit.cpp:148
void deployResults(double *vals, std::vector< int > *removable_ids)
Definition Circuit.cpp:491
double getResistance(std::string name)
Definition Circuit.cpp:90
double alphaBest
Best alpha scaling value.
Definition Circuit.h:86
bool checkCircuit(std::string substationId="")
Definition Circuit.cpp:905
std::vector< Element * > * voltageSources
Definition Circuit.h:68
@ ALPHA_VOLTAGE_LIMITS
The scaling alpha is applied (is not one] due to voltage limits.
Definition Circuit.h:102
@ ALPHA_NOT_APPLIED
The scaling alpha is not applied (is one)
Definition Circuit.h:98
@ ALPHA_CURRENT_LIMITS
The scaling alpha is applied (is not one) due to current limits.
Definition Circuit.h:100
@ ALPHA_NOT_CONVERGING
The Newton-Rhapson method has reached maximum iterations and no solution of circuit has been found wi...
Definition Circuit.h:104
Node * getNode(std::string name)
Definition Circuit.cpp:98
void descreaseLastId()
Definition Circuit.h:245
bool iscleaned
Definition Circuit.h:71
std::string & getCurrentsOfCircuitSource(std::string &currents)
List of currents of voltage sources as a string.
Definition Circuit.cpp:165
void eraseElement(Element *element)
Definition Circuit.cpp:834
void setId(int id)
Definition Element.cpp:133
ElementType getType()
Definition Element.cpp:119
double getResistance()
Definition Element.cpp:99
void setNegNode(Node *node)
Definition Element.cpp:130
Node * getNegNode()
Definition Element.cpp:115
double getCurrent()
Definition Element.cpp:85
void setPosNode(Node *node)
Definition Element.cpp:126
Node * getPosNode()
Definition Element.cpp:112
double getVoltage()
Definition Element.cpp:76
bool isEnabled()
Definition Element.cpp:148
double getPower()
Definition Element.cpp:105
Node * getTheOtherNode(Node *node)
Definition Element.cpp:138
ElementType
Definition Element.h:53
@ RESISTOR_traction_wire
Definition Element.h:54
@ CURRENT_SOURCE_traction_wire
Definition Element.h:55
@ VOLTAGE_SOURCE_traction_wire
Definition Element.h:56
@ ERROR_traction_wire
Definition Element.h:57
static bool gOverheadWireCurrentLimits
Definition MSGlobals.h:124
Definition Node.h:39
void setVoltage(double volt)
Definition Node.cpp:57
void setGround(bool isground)
Definition Node.cpp:73
void setNumMatrixCol(int num)
Definition Node.cpp:93
int getId()
Definition Node.cpp:77
Element * getAnOtherElement(Element *element)
Definition Node.cpp:109
bool isGround()
Definition Node.cpp:69
void addElement(Element *element)
Definition Node.cpp:44
void eraseElement(Element *element)
Definition Node.cpp:48
double getVoltage()
Definition Node.cpp:53
void setId(int id)
Definition Node.cpp:81
std::vector< Element * > * getElements()
Definition Node.cpp:101
bool isRemovable() const
Definition Node.h:68