Line data Source code
1 : /****************************************************************************/
2 : // Eclipse SUMO, Simulation of Urban MObility; see https://eclipse.dev/sumo
3 : // Copyright (C) 2001-2024 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 1828 : for (int i = 0; i < numofeqs - (int) voltageSources->size(); i++) {
299 1387 : vals[i] = 0;
300 : }
301 : J = A;
302 :
303 : int i = 0;
304 6963 : for (auto& node : *nodes) {
305 6522 : if (node->isGround() || node->isRemovable() || node->getNumMatrixRow() == -2) {
306 5135 : continue;
307 : }
308 1387 : 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 3740 : for (auto it_element = node->getElements()->begin(); it_element != node->getElements()->end(); it_element++) {
314 2353 : if ((*it_element)->getType() == Element::ElementType::CURRENT_SOURCE_traction_wire) {
315 : // if the element is current source
316 441 : if ((*it_element)->isEnabled()) {
317 : double diff_voltage;
318 441 : int PosNode_NumACol = (*it_element)->getPosNode()->getNumMatrixCol();
319 441 : int NegNode_NumACol = (*it_element)->getNegNode()->getNumMatrixCol();
320 : // compute voltage on current source
321 441 : if (PosNode_NumACol == -1) {
322 : // if the positive node is the ground => U = 0 - phi(NegNode)
323 0 : diff_voltage = -x[NegNode_NumACol];
324 441 : } else if (NegNode_NumACol == -1) {
325 : // if the negative node is the ground => U = phi(PosNode) - 0
326 441 : 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 441 : 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 441 : vals[i] -= alpha * (*it_element)->getPowerWanted() / diff_voltage;
335 441 : (*it_element)->setCurrent(-alpha * (*it_element)->getPowerWanted() / diff_voltage);
336 441 : 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 441 : J(i, PosNode_NumACol) -= alpha * (*it_element)->getPowerWanted() / diff_voltage / diff_voltage;
339 : }
340 441 : 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 1387 : 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 1828 : for (i = 0; i < numofeqs - (int)voltageSources->size(); i++) {
374 1387 : currentSumActual -= vals[i];
375 : }
376 : // RICE_TODO @20210409 This epsilon should be specified somewhere as a constant. Or should be a parameter.
377 441 : 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 261 : } 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 261 : dx = -J.colPivHouseholderQr().solve(A * x - b);
413 261 : x = x + dx;
414 261 : ++iterNR;
415 261 : }
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 8 : Circuit::Circuit(double currentLimit) {
591 8 : nodes = new std::vector<Node*>(0);
592 8 : elements = new std::vector<Element*>(0);
593 8 : voltageSources = new std::vector<Element*>(0);
594 8 : lastId = 0;
595 8 : iscleaned = true;
596 8 : circuitCurrentLimit = currentLimit;
597 8 : }
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 : break;
747 : case Element::ElementType::ERROR_traction_wire:
748 : return false;
749 : break;
750 : }
751 : }
752 : return true;
753 : }
754 : #endif
755 :
756 : /**
757 : * Select removable nodes, i.e. nodes that are NOT the ground of the circuit
758 : * and that have exactly two resistor elements connected. Ids of those
759 : * removable nodes are added into the internal vector `removable_ids`.
760 : */
761 180 : void Circuit::detectRemovableNodes(std::vector<int>* removable_ids) {
762 : // loop over all nodes in the circuit
763 2793 : for (std::vector<Node*>::iterator it = nodes->begin(); it != nodes->end(); it++) {
764 : // if the node is connected to two elements and is not the ground
765 2613 : if ((*it)->getElements()->size() == 2 && !(*it)->isGround()) {
766 : // set such node by default as removable. But check if the two elements are both resistors
767 1691 : (*it)->setRemovability(true);
768 5031 : for (std::vector<Element*>::iterator it2 = (*it)->getElements()->begin(); it2 != (*it)->getElements()->end(); it2++) {
769 3361 : if ((*it2)->getType() != Element::ElementType::RESISTOR_traction_wire) {
770 21 : (*it)->setRemovability(false);
771 : break;
772 : }
773 : }
774 1691 : if ((*it)->isRemovable()) {
775 : //if the node is removable add pointer into the vector of removable nodes
776 1670 : removable_ids->push_back((*it)->getId());
777 : }
778 : } else {
779 922 : (*it)->setRemovability(false);
780 : }
781 : }
782 : // sort the vector of removable ids
783 180 : std::sort(removable_ids->begin(), removable_ids->end(), std::less<int>());
784 180 : return;
785 : }
786 :
787 424 : Element* Circuit::addElement(std::string name, double value, Node* pNode, Node* nNode, Element::ElementType et) {
788 : // RICE_CHECK: This seems to be a bit of work in progress, is it final?
789 : // if ((et == Element::ElementType::RESISTOR_traction_wire && value <= 0) || et == Element::ElementType::ERROR_traction_wire) {
790 424 : if (et == Element::ElementType::RESISTOR_traction_wire && value <= 1e-6) {
791 : //due to numeric problems
792 : // RICE_TODO @20210409 This epsilon should be specified somewhere as a constant. Or should be a parameter.
793 0 : if (value > -1e-6) {
794 : value = 1e-6;
795 0 : WRITE_WARNING(TL("Trying to add resistor element into the overhead wire circuit with resistance < 1e-6. "))
796 : } else {
797 0 : WRITE_ERROR(TL("Trying to add resistor element into the overhead wire circuit with resistance < 0. "))
798 0 : return nullptr;
799 : }
800 : }
801 :
802 424 : Element* e = getElement(name);
803 :
804 424 : if (e != nullptr) {
805 : //WRITE_ERRORF(TL("The element: '%' already exists."), name);
806 0 : std::cout << "The element: '" + name + "' already exists.";
807 0 : return nullptr;
808 : }
809 :
810 848 : e = new Element(name, et, value);
811 424 : if (e->getType() == Element::ElementType::VOLTAGE_SOURCE_traction_wire) {
812 11 : e->setId(lastId);
813 11 : lastId++;
814 11 : circuit_lock.lock();
815 11 : this->voltageSources->push_back(e);
816 : circuit_lock.unlock();
817 : } else {
818 413 : circuit_lock.lock();
819 413 : this->elements->push_back(e);
820 : circuit_lock.unlock();
821 : }
822 :
823 424 : e->setPosNode(pNode);
824 424 : e->setNegNode(nNode);
825 :
826 424 : pNode->addElement(e);
827 424 : nNode->addElement(e);
828 : return e;
829 : }
830 :
831 402 : void Circuit::eraseElement(Element* element) {
832 402 : element->getPosNode()->eraseElement(element);
833 402 : element->getNegNode()->eraseElement(element);
834 402 : circuit_lock.lock();
835 402 : this->elements->erase(std::remove(this->elements->begin(), this->elements->end(), element), this->elements->end());
836 : circuit_lock.unlock();
837 402 : }
838 :
839 15 : void Circuit::replaceAndDeleteNode(Node* unusedNode, Node* newNode) {
840 : //replace element node if it is unusedNode
841 27 : for (auto& voltageSource : *voltageSources) {
842 12 : if (voltageSource->getNegNode() == unusedNode) {
843 0 : voltageSource->setNegNode(newNode);
844 0 : newNode->eraseElement(voltageSource);
845 0 : newNode->addElement(voltageSource);
846 : }
847 12 : if (voltageSource->getPosNode() == unusedNode) {
848 0 : voltageSource->setPosNode(newNode);
849 0 : newNode->eraseElement(voltageSource);
850 0 : newNode->addElement(voltageSource);
851 : }
852 : }
853 60 : for (auto& element : *elements) {
854 45 : if (element->getNegNode() == unusedNode) {
855 0 : element->setNegNode(newNode);
856 0 : newNode->eraseElement(element);
857 0 : newNode->addElement(element);
858 : }
859 45 : if (element->getPosNode() == unusedNode) {
860 15 : element->setPosNode(newNode);
861 15 : newNode->eraseElement(element);
862 15 : newNode->addElement(element);
863 : }
864 : }
865 :
866 : //erase unusedNode from nodes vector
867 15 : this->eraseNode(unusedNode);
868 :
869 : //modify id of other elements and nodes
870 15 : int modLastId = this->getLastId() - 1;
871 15 : if (unusedNode->getId() != modLastId) {
872 15 : Node* node_last = this->getNode(modLastId);
873 15 : if (node_last != nullptr) {
874 15 : node_last->setId(unusedNode->getId());
875 : } else {
876 0 : Element* elem_last = this->getVoltageSource(modLastId);
877 0 : if (elem_last != nullptr) {
878 0 : elem_last->setId(unusedNode->getId());
879 : } else {
880 0 : WRITE_ERROR(TL("The element or node with the last Id was not found in the circuit!"));
881 : }
882 : }
883 : }
884 :
885 : this->decreaseLastId();
886 15 : delete unusedNode;
887 15 : }
888 :
889 0 : void Circuit::cleanUpSP() {
890 0 : for (std::vector<Element*>::iterator it = elements->begin(); it != elements->end(); it++) {
891 0 : if ((*it)->getType() != Element::ElementType::RESISTOR_traction_wire) {
892 0 : (*it)->setEnabled(true);
893 : }
894 : }
895 :
896 0 : for (std::vector<Element*>::iterator it = voltageSources->begin(); it != voltageSources->end(); it++) {
897 0 : (*it)->setEnabled(true);
898 : }
899 0 : this->iscleaned = true;
900 0 : }
901 :
902 8 : bool Circuit::checkCircuit(std::string substationId) {
903 : // check empty nodes
904 77 : for (std::vector<Node*>::iterator it = nodes->begin(); it != nodes->end(); it++) {
905 69 : if ((*it)->getNumOfElements() < 2) {
906 : //cout << "WARNING: Node [" << (*it)->getName() << "] is connected to less than two elements, please enter other elements.\n";
907 24 : if ((*it)->getNumOfElements() < 1) {
908 : return false;
909 : }
910 : }
911 : }
912 : // check voltage sources
913 19 : for (std::vector<Element*>::iterator it = voltageSources->begin(); it != voltageSources->end(); it++) {
914 11 : if ((*it)->getPosNode() == nullptr || (*it)->getNegNode() == nullptr) {
915 : //cout << "ERROR: Voltage Source [" << (*it)->getName() << "] is connected to less than two nodes, please enter the other end.\n";
916 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);
917 : return false;
918 : }
919 : }
920 : // check other elements
921 61 : for (std::vector<Element*>::iterator it = elements->begin(); it != elements->end(); it++) {
922 53 : if ((*it)->getPosNode() == nullptr || (*it)->getNegNode() == nullptr) {
923 : //cout << "ERROR: Element [" << (*it)->getName() << "] is connected to less than two nodes, please enter the other end.\n";
924 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);
925 : return false;
926 : }
927 : }
928 :
929 : // check connectivity
930 8 : int num = (int)nodes->size() + getNumVoltageSources() - 1;
931 8 : bool* nodesVisited = new bool[num];
932 80 : for (int i = 0; i < num; i++) {
933 72 : nodesVisited[i] = false;
934 : }
935 : // TODO: Probably unused
936 : // int id = -1;
937 8 : if (!getNode(-1)->isGround()) {
938 : //cout << "ERROR: Node id -1 is not the ground \n";
939 0 : WRITE_ERRORF(TL("Circuit Node with id '-1' is not the grounded, please adjust the definition of the section (with substation '%')."), substationId);
940 : }
941 8 : std::vector<Node*>* queue = new std::vector<Node*>(0);
942 8 : Node* node = nullptr;
943 8 : Node* neigboringNode = nullptr;
944 : //start with (voltageSources->front()->getPosNode())
945 8 : nodesVisited[voltageSources->front()->getId()] = 1;
946 8 : node = voltageSources->front()->getPosNode();
947 8 : queue->push_back(node);
948 :
949 122 : while (!queue->empty()) {
950 114 : node = queue->back();
951 : queue->pop_back();
952 114 : if (!nodesVisited[node->getId()]) {
953 61 : nodesVisited[node->getId()] = true;
954 178 : for (auto it = node->getElements()->begin(); it != node->getElements()->end(); it++) {
955 117 : neigboringNode = (*it)->getTheOtherNode(node);
956 117 : if (!neigboringNode->isGround()) {
957 106 : queue->push_back(neigboringNode);
958 11 : } else if ((*it)->getType() == Element::ElementType::VOLTAGE_SOURCE_traction_wire) {
959 : /// there used to be == 1 which was probably a typo ... check!
960 11 : nodesVisited[(*it)->getId()] = 1;
961 0 : } else if ((*it)->getType() == Element::ElementType::RESISTOR_traction_wire) {
962 : //cout << "ERROR: The resistor type connects the ground \n";
963 0 : WRITE_ERRORF(TL("A Circuit Resistor Element connects the ground, please adjust the definition of the section (with substation '%')."), substationId);
964 : }
965 : }
966 : }
967 : }
968 :
969 80 : for (int i = 0; i < num; i++) {
970 72 : if (nodesVisited[i] == 0) {
971 : //cout << "ERROR: Node or voltage source with id " << (i) << " has been not visited during checking of the circuit => Disconnectivity of the circuit. \n";
972 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);
973 : }
974 : }
975 :
976 8 : return true;
977 : }
978 :
979 196 : int Circuit::getNumVoltageSources() {
980 196 : return (int) voltageSources->size();
981 : }
|
