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