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