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