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