39 #pragma warning(disable: 4127)
41 #include "Eigen/Dense"
42 #include "Eigen/Sparse"
43 #include "Eigen/Geometry"
57 WRITE_ERROR(
"The node: '" + name +
"' already exists.");
61 if (
nodes->size() == 0) {
70 this->
nodes->push_back(tNode);
77 this->
nodes->erase(std::remove(this->
nodes->begin(), this->nodes->end(), node), this->nodes->end());
83 if (tElement ==
nullptr) {
91 if (tElement ==
nullptr) {
93 if (node !=
nullptr) {
105 if (tElement ==
nullptr) {
113 if (node->getName() == name) {
122 if (node->getId() ==
id) {
131 if (el->getName() == name) {
136 if (voltageSource->getName() == name) {
137 return voltageSource;
145 if (el->getId() ==
id) {
154 if (voltageSource->getId() ==
id) {
155 return voltageSource;
164 power += voltageSource->getPower();
172 current += voltageSource->getCurrent();
182 currents +=
toString(voltageSource->getCurrent(), 4) +
" ";
184 if (!currents.empty()) {
192 std::vector<Element*>* vsources =
new std::vector<Element*>(0);
194 if (el->getType() == Element::ElementType::CURRENT_SOURCE_traction_wire) {
196 vsources->push_back(el);
211 void Circuit::removeColumn(Eigen::MatrixXd& matrix,
int colToRemove) {
212 const int numRows = (int)matrix.rows();
213 const int numCols = (int)matrix.cols() - 1;
215 if (colToRemove < numCols) {
216 matrix.block(0, colToRemove, numRows, numCols - colToRemove) = matrix.rightCols(numCols - colToRemove);
219 matrix.conservativeResize(numRows, numCols);
222 bool Circuit::solveEquationsNRmethod(
double* eqn,
double* vals, std::vector<int>* removable_ids) {
225 int numofeqs = numofcolumn - (int)removable_ids->size();
228 Eigen::MatrixXd A = Eigen::Map < Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> >(eqn, numofeqs, numofcolumn);
233 for (std::vector<int>::reverse_iterator it = removable_ids->rbegin(); it != removable_ids->rend(); ++it) {
234 id = (*it >= 0 ? *it : -(*it));
243 Node* tNode =
nullptr;
244 for (
int i = 0; i < numofcolumn; i++) {
246 if (tNode !=
nullptr) {
252 WRITE_ERROR(
"Index of renumbered node exceeded the reduced number of equations.");
261 if (tElem !=
nullptr) {
263 WRITE_ERROR(
"Index of renumbered element exceeded the reduced number of equations.");
270 WRITE_ERROR(
"Structural error in reduced circuit matrix.");
274 Eigen::Map<Eigen::VectorXd> b(vals, numofeqs);
275 Eigen::VectorXd x = A.colPivHouseholderQr().solve(b);
278 Eigen::MatrixXd J = A;
281 int max_iter_of_NR = 10;
291 std::vector<double> alpha_notSolution;
293 double alpha_res = 1e-2;
295 double currentSumActual = 0.0;
297 Eigen::VectorXd x_best = x;
298 bool x_best_exist =
true;
300 if (x.maxCoeff() > 10e6 || x.minCoeff() < -10e6) {
301 WRITE_ERROR(
"Initial solution x used during solving DC circuit is out of bounds.\n");
314 for (
int i = 0; i < numofeqs - (int)
voltageSources->size(); i++) {
320 for (
auto& node : *
nodes) {
321 if (node->isGround() || node->isRemovable() || node->getNumMatrixRow() == -2) {
324 if (node->getNumMatrixRow() != i) {
325 WRITE_ERROR(
"wrongly assigned row of matrix A during solving the circuit");
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) {
332 if ((*it_element)->isEnabled()) {
334 int PosNode_NumACol = (*it_element)->getPosNode()->getNumMatrixCol();
335 int NegNode_NumACol = (*it_element)->getNegNode()->getNumMatrixCol();
337 if (PosNode_NumACol == -1) {
339 diff_voltage = -x[NegNode_NumACol];
340 }
else if (NegNode_NumACol == -1) {
342 diff_voltage = x[PosNode_NumACol];
345 diff_voltage = (x[PosNode_NumACol] - x[NegNode_NumACol]);
348 if ((*it_element)->getPosNode() == node) {
350 vals[i] -= alpha * (*it_element)->getPowerWanted() / diff_voltage;
351 (*it_element)->setCurrent(-alpha * (*it_element)->getPowerWanted() / diff_voltage);
352 if (PosNode_NumACol != -1) {
354 J(i, PosNode_NumACol) -= alpha * (*it_element)->getPowerWanted() / diff_voltage / diff_voltage;
356 if (NegNode_NumACol != -1) {
358 J(i, NegNode_NumACol) += alpha * (*it_element)->getPowerWanted() / diff_voltage / diff_voltage;
362 vals[i] += alpha * (*it_element)->getPowerWanted() / diff_voltage;
367 WRITE_WARNING(
"The negative node of current source is not the groud.")
368 if (PosNode_NumACol != -1) {
370 J(i, PosNode_NumACol) += alpha * (*it_element)->getPowerWanted() / diff_voltage / diff_voltage;
372 if (NegNode_NumACol != -1) {
374 J(i, NegNode_NumACol) -= alpha * (*it_element)->getPowerWanted() / diff_voltage / diff_voltage;
388 currentSumActual = 0;
390 currentSumActual -= vals[i];
393 if ((A * x - b).norm() < 1e-6) {
397 alpha_notSolution.push_back(alpha);
407 alpha_notSolution.push_back(alpha);
418 }
else if (iterNR == max_iter_of_NR) {
420 alpha_notSolution.push_back(alpha);
428 dx = -J.colPivHouseholderQr().solve(A * x - b);
433 if (alpha_notSolution.empty()) {
438 if ((alpha_notSolution.back() -
alphaBest) < alpha_res) {
439 max_iter_of_NR = 2 * max_iter_of_NR;
443 alpha_res = alpha_res / 10;
445 if (alpha_res < 5e-5) {
448 alpha = alpha_notSolution.back();
449 alpha_notSolution.pop_back();
457 for (
int i = 0; i < numofeqs; i++) {
464 for (
auto& node : *
nodes) {
465 if (node->isGround() || node->isRemovable() || node->getNumMatrixRow() == -2) {
468 if (node->getNumMatrixRow() != i) {
469 WRITE_ERROR(
"wrongly assigned row of matrix A during solving the circuit");
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()) {
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];
482 diff_voltage = (x_best[PosNode_NumACol] - x_best[NegNode_NumACol]);
485 if ((*it_element)->getPosNode() == node) {
486 (*it_element)->setCurrent(-
alphaBest * (*it_element)->getPowerWanted() / diff_voltage);
492 WRITE_WARNING(
"The negative node of current source is not the groud.")
508 int numofeqs = numofcolumn - (int)removable_ids->size();
513 Node* tNode =
nullptr;
514 for (
int i = 0; i < numofcolumn; i++) {
516 if (tNode !=
nullptr)
521 WRITE_ERROR(
"Results deployment during circuit evaluation was unsuccessfull.");
529 if (tElem !=
nullptr) {
531 WRITE_ERROR(
"Results deployment during circuit evaluation was unsuccessfull.");
539 WRITE_ERROR(
"Results deployment during circuit evaluation was unsuccessfull.");
544 Node* nextNONremovableNode1 =
nullptr;
545 Node* nextNONremovableNode2 =
nullptr;
548 if (!node->isRemovable()) {
551 if (node->getElements()->size() != 2) {
555 el1 = node->getElements()->front();
556 el2 = node->getElements()->back();
575 y = ((1 - x) * nextNONremovableNode1->
getVoltage()) + (x * nextNONremovableNode2->
getVoltage());
577 node->setRemovability(
false);
582 double currentSum = 0;
585 if (el != voltageSource) {
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.");
593 voltageSource->setCurrent(currentSum);
598 nodes =
new std::vector<Node*>(0);
599 elements =
new std::vector<Element*>(0);
607 nodes =
new std::vector<Node*>(0);
608 elements =
new std::vector<Element*>(0);
616 bool Circuit::_solveNRmethod() {
617 double* eqn =
nullptr;
618 double* vals =
nullptr;
619 std::vector<int> removable_ids;
622 createEquationsNRmethod(eqn, vals, &removable_ids);
623 if (!solveEquationsNRmethod(eqn, vals, &removable_ids)) {
634 bool Circuit::solve() {
638 return this->_solveNRmethod();
641 bool Circuit::createEquationsNRmethod(
double*& eqs,
double*& vals, std::vector<int>* removable_ids) {
649 int m = n - (int)(removable_ids->size() +
voltageSources->size());
652 eqs =
new double[m * n];
653 vals =
new double[m];
655 for (
int i = 0; i < m; i++) {
657 for (
int j = 0; j < n; j++) {
664 for (std::vector<Node*>::iterator it =
nodes->begin(); it !=
nodes->end(); it++) {
665 if ((*it)->isGround() || (*it)->isRemovable()) {
667 (*it)->setNumMatrixRow(-1);
672 bool noVoltageSource = createEquationNRmethod((*it), (eqs + n * i), vals[i], removable_ids);
674 if (noVoltageSource) {
675 (*it)->setNumMatrixRow(i);
678 (*it)->setNumMatrixRow(-2);
680 for (
int j = 0; j < n; j++) {
687 std::sort(removable_ids->begin(), removable_ids->end(), std::less<int>());
692 createEquation((*it), (eqs + n * i), vals[i]);
699 bool Circuit::createEquation(
Element* vsource,
double* eqn,
double& val) {
714 bool Circuit::createEquationNRmethod(
Node* node,
double* eqn,
double& val, std::vector<int>* removable_ids) {
716 for (std::vector<Element*>::iterator it = node->
getElements()->begin(); it != node->
getElements()->end(); it++) {
718 switch ((*it)->getType()) {
719 case Element::ElementType::RESISTOR_traction_wire:
720 if ((*it)->isEnabled()) {
721 x = (*it)->getResistance();
723 Node* nextNONremovableNode = (*it)->getTheOtherNode(node);
724 Element* nextSerialResistor = *it;
726 nextSerialResistor = nextNONremovableNode->
getAnOtherElement(nextSerialResistor);
728 nextNONremovableNode = nextSerialResistor->
getTheOtherNode(nextNONremovableNode);
732 eqn[node->
getId()] += x;
734 if (!nextNONremovableNode->
isGround()) {
735 eqn[nextNONremovableNode->
getId()] -= x;
739 case Element::ElementType::CURRENT_SOURCE_traction_wire:
740 if ((*it)->isEnabled()) {
742 if ((*it)->getPosNode() == node) {
743 x = -(*it)->getPowerWanted() /
voltageSources->front()->getVoltage();
745 x = (*it)->getPowerWanted() /
voltageSources->front()->getVoltage();
752 case Element::ElementType::VOLTAGE_SOURCE_traction_wire:
753 if ((*it)->getPosNode() == node) {
758 eqn[(*it)->getId()] += x;
760 removable_ids->push_back((*it)->getId());
763 case Element::ElementType::ERROR_traction_wire:
779 for (std::vector<Node*>::iterator it =
nodes->begin(); it !=
nodes->end(); it++) {
781 if ((*it)->getElements()->size() == 2 && !(*it)->isGround()) {
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);
790 if ((*it)->isRemovable()) {
792 removable_ids->push_back((*it)->getId());
795 (*it)->setRemovability(
false);
799 std::sort(removable_ids->begin(), removable_ids->end(), std::less<int>());
806 if (et == Element::ElementType::RESISTOR_traction_wire && value <= 1e-6) {
811 WRITE_WARNING(
"Trying to add resistor element into the overhead wire circuit with resistance < 1e-6. ")
813 WRITE_ERROR(
"Trying to add resistor element into the overhead wire circuit with resistance < 0. ")
822 std::cout <<
"The element: '" + name +
"' already exists.";
826 e =
new Element(name, et, value);
827 if (e->
getType() == Element::ElementType::VOLTAGE_SOURCE_traction_wire) {
851 this->
elements->erase(std::remove(this->
elements->begin(), this->elements->end(), element), this->elements->end());
858 if (voltageSource->getNegNode() == unusedNode) {
859 voltageSource->setNegNode(newNode);
863 if (voltageSource->getPosNode() == unusedNode) {
864 voltageSource->setPosNode(newNode);
870 if (element->getNegNode() == unusedNode) {
871 element->setNegNode(newNode);
875 if (element->getPosNode() == unusedNode) {
876 element->setPosNode(newNode);
887 if (unusedNode->
getId() != modLastId) {
889 if (node_last !=
nullptr) {
893 if (elem_last !=
nullptr) {
896 WRITE_ERROR(
"The element or node with the last Id was not found in the circuit!");
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);
913 (*it)->setEnabled(
true);
920 for (std::vector<Node*>::iterator it =
nodes->begin(); it !=
nodes->end(); it++) {
921 if ((*it)->getNumOfElements() < 2) {
923 if ((*it)->getNumOfElements() < 1) {
930 if ((*it)->getPosNode() ==
nullptr || (*it)->getNegNode() ==
nullptr) {
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 +
"').");
937 for (std::vector<Element*>::iterator it =
elements->begin(); it !=
elements->end(); it++) {
938 if ((*it)->getPosNode() ==
nullptr || (*it)->getNegNode() ==
nullptr) {
940 WRITE_ERROR(
"Circuit Element '" + (*it)->getName() +
"' is connected to less than two nodes, please adjust the definition of the section (with substation '" + substationId +
"').");
947 bool* nodesVisited =
new bool[num];
948 for (
int i = 0; i < num; i++) {
949 nodesVisited[i] =
false;
953 if (!
getNode(-1)->isGround()) {
955 WRITE_ERROR(
"Circuit Node with id '-1' is not the grounded, please adjust the definition of the section (with substation '" + substationId +
"').");
957 std::vector<Node*>* queue =
new std::vector<Node*>(0);
958 Node* node =
nullptr;
959 Node* neigboringNode =
nullptr;
963 queue->push_back(node);
965 while (!queue->empty()) {
966 node = queue->back();
968 if (!nodesVisited[node->
getId()]) {
969 nodesVisited[node->
getId()] =
true;
971 neigboringNode = (*it)->getTheOtherNode(node);
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) {
979 WRITE_ERROR(
"A Circuit Resistor Element connects the ground, please adjust the definition of the section (with substation '" + substationId +
"').");
985 for (
int i = 0; i < num; i++) {
986 if (nodesVisited[i] == 0) {
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 +
"').");
static std::mutex circuit_lock
#define WRITE_WARNING(msg)
std::string toString(const T &t, std::streamsize accuracy=gPrecision)
std::vector< Node * > * nodes
std::vector< Element * > * getCurrentSources()
Node * addNode(std::string name)
double getVoltage(std::string name)
int getNumVoltageSources()
Element * addElement(std::string name, double value, Node *pNode, Node *nNode, Element::ElementType et)
void eraseNode(Node *node)
double getCurrent(std::string name)
Element * getElement(std::string name)
double circuitCurrentLimit
The electric current limit of the voltage sources.
double getTotalCurrentOfCircuitSources()
The sum of voltage source currents in the circuit.
void detectRemovableNodes(std::vector< int > *removable_ids)
double getCurrentLimit()
@ brief Get the electric current limit of this circuit.
std::vector< Element * > * elements
void replaceAndDeleteNode(Node *unusedNode, Node *newNode)
Element * getVoltageSource(int id)
double getTotalPowerOfCircuitSources()
The sum of voltage source powers in the circuit.
void deployResults(double *vals, std::vector< int > *removable_ids)
double getResistance(std::string name)
double alphaBest
Best alpha scaling value.
bool checkCircuit(std::string substationId="")
std::vector< Element * > * voltageSources
@ ALPHA_VOLTAGE_LIMITS
The scaling alpha is applied (is not one] due to voltage limits.
@ ALPHA_NOT_APPLIED
The scaling alpha is not applied (is one)
@ ALPHA_CURRENT_LIMITS
The scaling alpha is applied (is not one) due to current limits.
@ ALPHA_NOT_CONVERGING
The Newton-Rhapson method has reached maximum iterations and no solution of circuit has been found wi...
Node * getNode(std::string name)
std::string & getCurrentsOfCircuitSource(std::string ¤ts)
List of currents of voltage sources as a string.
void eraseElement(Element *element)
void setNegNode(Node *node)
void setPosNode(Node *node)
Node * getTheOtherNode(Node *node)
static bool gOverheadWireCurrentLimits
void setVoltage(double volt)
std::vector< Element * > * getElements()
void setGround(bool isground)
void setNumMatrixCol(int num)
Element * getAnOtherElement(Element *element)
void addElement(Element *element)
void eraseElement(Element *element)