48 if (
nodes->size() == 0) {
57 this->
nodes->push_back(tNode);
64 this->
nodes->erase(std::remove(this->
nodes->begin(), this->nodes->end(), node), this->nodes->end());
70 if (tElement ==
nullptr) {
78 if (tElement ==
nullptr) {
80 if (node !=
nullptr) {
92 if (tElement ==
nullptr) {
100 if (node->getName() == name) {
109 if (node->getId() == id) {
118 if (el->getName() == name) {
123 if (voltageSource->getName() == name) {
124 return voltageSource;
132 if (el->getId() == id) {
141 if (voltageSource->getId() == id) {
142 return voltageSource;
159 current += voltageSource->getCurrent();
169 currents +=
toString(voltageSource->getCurrent(), 4) +
" ";
171 if (!currents.empty()) {
179 std::vector<Element*>* vsources =
new std::vector<Element*>(0);
183 vsources->push_back(el);
198void Circuit::removeColumn(Eigen::MatrixXd& matrix,
int colToRemove) {
199 const int numRows = (int)matrix.rows();
200 const int numCols = (int)matrix.cols() - 1;
202 if (colToRemove < numCols) {
203 matrix.block(0, colToRemove, numRows, numCols - colToRemove) = matrix.rightCols(numCols - colToRemove);
206 matrix.conservativeResize(numRows, numCols);
209bool Circuit::solveEquationsNRmethod(
double* eqn,
double* vals, std::vector<int>* removable_ids) {
212 int numofeqs = numofcolumn - (int)removable_ids->size();
215 Eigen::MatrixXd A = Eigen::Map < Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> >(eqn, numofeqs, numofcolumn);
220 for (std::vector<int>::reverse_iterator it = removable_ids->rbegin(); it != removable_ids->rend(); ++it) {
221 id = (*it >= 0 ? *it : -(*it));
230 Node* tNode =
nullptr;
231 for (
int i = 0; i < numofcolumn; i++) {
233 if (tNode !=
nullptr) {
239 WRITE_ERROR(
TL(
"Index of renumbered node exceeded the reduced number of equations."));
248 if (tElem !=
nullptr) {
250 WRITE_ERROR(
TL(
"Index of renumbered element exceeded the reduced number of equations."));
257 WRITE_ERROR(
TL(
"Structural error in reduced circuit matrix."));
261 Eigen::Map<Eigen::VectorXd> b(vals, numofeqs);
262 Eigen::VectorXd x = A.colPivHouseholderQr().solve(b);
265 Eigen::MatrixXd J = A;
268 int max_iter_of_NR = 10;
278 std::vector<double> alpha_notSolution;
280 double alpha_res = 1e-2;
282 double currentSumActual = 0.0;
284 Eigen::VectorXd x_best = x;
285 bool x_best_exist =
true;
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"));
301 for (
int i = 0; i < numofeqs - (int)
voltageSources->size(); i++) {
307 for (
auto& node : *
nodes) {
308 if (node->isGround() || node->isRemovable() || node->getNumMatrixRow() == -2) {
311 if (node->getNumMatrixRow() != i) {
312 WRITE_ERROR(
TL(
"wrongly assigned row of matrix A during solving the circuit"));
316 for (
auto it_element = node->getElements()->begin(); it_element != node->getElements()->end(); it_element++) {
319 if ((*it_element)->isEnabled()) {
321 int PosNode_NumACol = (*it_element)->getPosNode()->getNumMatrixCol();
322 int NegNode_NumACol = (*it_element)->getNegNode()->getNumMatrixCol();
324 if (PosNode_NumACol == -1) {
326 diff_voltage = -x[NegNode_NumACol];
327 }
else if (NegNode_NumACol == -1) {
329 diff_voltage = x[PosNode_NumACol];
332 diff_voltage = (x[PosNode_NumACol] - x[NegNode_NumACol]);
335 if ((*it_element)->getPosNode() == node) {
337 vals[i] -= alpha * (*it_element)->getPowerWanted() / diff_voltage;
338 (*it_element)->setCurrent(-alpha * (*it_element)->getPowerWanted() / diff_voltage);
339 if (PosNode_NumACol != -1) {
341 J(i, PosNode_NumACol) -= alpha * (*it_element)->getPowerWanted() / diff_voltage / diff_voltage;
343 if (NegNode_NumACol != -1) {
345 J(i, NegNode_NumACol) += alpha * (*it_element)->getPowerWanted() / diff_voltage / diff_voltage;
349 vals[i] += alpha * (*it_element)->getPowerWanted() / diff_voltage;
354 WRITE_WARNING(
TL(
"The negative node of current source is not the groud."))
355 if (PosNode_NumACol != -1) {
357 J(i, PosNode_NumACol) += alpha * (*it_element)->getPowerWanted() / diff_voltage / diff_voltage;
359 if (NegNode_NumACol != -1) {
361 J(i, NegNode_NumACol) -= alpha * (*it_element)->getPowerWanted() / diff_voltage / diff_voltage;
375 currentSumActual = 0;
377 currentSumActual -= vals[i];
380 if ((A * x - b).norm() < 1e-6) {
384 alpha_notSolution.push_back(alpha);
394 alpha_notSolution.push_back(alpha);
405 }
else if (iterNR == max_iter_of_NR) {
407 alpha_notSolution.push_back(alpha);
415 dx = -J.colPivHouseholderQr().solve(A * x - b);
420 if (alpha_notSolution.empty()) {
425 if ((alpha_notSolution.back() -
alphaBest) < alpha_res) {
426 max_iter_of_NR = 2 * max_iter_of_NR;
430 alpha_res = alpha_res / 10;
432 if (alpha_res < 5e-5) {
435 alpha = alpha_notSolution.back();
436 alpha_notSolution.pop_back();
444 for (
int i = 0; i < numofeqs; i++) {
451 for (
auto& node : *
nodes) {
452 if (node->isGround() || node->isRemovable() || node->getNumMatrixRow() == -2) {
455 if (node->getNumMatrixRow() != i) {
456 WRITE_ERROR(
TL(
"wrongly assigned row of matrix A during solving the circuit"));
458 for (
auto it_element = node->getElements()->begin(); it_element != node->getElements()->end(); it_element++) {
460 if ((*it_element)->isEnabled()) {
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];
469 diff_voltage = (x_best[PosNode_NumACol] - x_best[NegNode_NumACol]);
472 if ((*it_element)->getPosNode() == node) {
473 (*it_element)->setCurrent(-
alphaBest * (*it_element)->getPowerWanted() / diff_voltage);
479 WRITE_WARNING(
TL(
"The negative node of current source is not the groud."))
495 int numofeqs = numofcolumn - (int)removable_ids->size();
500 Node* tNode =
nullptr;
501 for (
int i = 0; i < numofcolumn; i++) {
503 if (tNode !=
nullptr)
508 WRITE_ERROR(
TL(
"Results deployment during circuit evaluation was unsuccessful."));
516 if (tElem !=
nullptr) {
518 WRITE_ERROR(
TL(
"Results deployment during circuit evaluation was unsuccessful."));
526 WRITE_ERROR(
TL(
"Results deployment during circuit evaluation was unsuccessful."));
531 Node* nextNONremovableNode1 =
nullptr;
532 Node* nextNONremovableNode2 =
nullptr;
535 if (!node->isRemovable()) {
538 if (node->getElements()->size() != 2) {
542 el1 = node->getElements()->front();
543 el2 = node->getElements()->back();
562 y = ((1 - x) * nextNONremovableNode1->
getVoltage()) + (x * nextNONremovableNode2->
getVoltage());
564 node->setRemovability(
false);
569 double currentSum = 0;
570 for (
Element*
const el : *voltageSource->getPosNode()->getElements()) {
572 if (el != voltageSource) {
574 currentSum += (voltageSource->getPosNode()->getVoltage() - el->getTheOtherNode(voltageSource->getPosNode())->
getVoltage()) / el->getResistance();
576 WRITE_WARNING(
TL(
"Cannot assign unambigous electric current value to two voltage sources connected in parallel at the same node."));
580 voltageSource->setCurrent(currentSum);
585 nodes =
new std::vector<Node*>(0);
586 elements =
new std::vector<Element*>(0);
594 nodes =
new std::vector<Node*>(0);
595 elements =
new std::vector<Element*>(0);
603bool Circuit::_solveNRmethod() {
604 double* eqn =
nullptr;
605 double* vals =
nullptr;
606 std::vector<int> removable_ids;
609 createEquationsNRmethod(eqn, vals, &removable_ids);
610 if (!solveEquationsNRmethod(eqn, vals, &removable_ids)) {
621bool Circuit::solve() {
625 return this->_solveNRmethod();
628bool Circuit::createEquationsNRmethod(
double*& eqs,
double*& vals, std::vector<int>* removable_ids) {
636 int m = n - (int)(removable_ids->size() +
voltageSources->size());
639 eqs =
new double[m * n];
640 vals =
new double[m];
642 for (
int i = 0; i < m; i++) {
644 for (
int j = 0; j < n; j++) {
651 for (std::vector<Node*>::iterator it =
nodes->begin(); it !=
nodes->end(); it++) {
652 if ((*it)->isGround() || (*it)->isRemovable()) {
654 (*it)->setNumMatrixRow(-1);
659 bool noVoltageSource = createEquationNRmethod((*it), (eqs + n * i), vals[i], removable_ids);
661 if (noVoltageSource) {
662 (*it)->setNumMatrixRow(i);
665 (*it)->setNumMatrixRow(-2);
667 for (
int j = 0; j < n; j++) {
674 std::sort(removable_ids->begin(), removable_ids->end(), std::less<int>());
679 createEquation((*it), (eqs + n * i), vals[i]);
686bool Circuit::createEquation(
Element* vsource,
double* eqn,
double& val) {
701bool Circuit::createEquationNRmethod(
Node* node,
double* eqn,
double& val, std::vector<int>* removable_ids) {
703 for (std::vector<Element*>::iterator it = node->
getElements()->begin(); it != node->
getElements()->end(); it++) {
705 switch ((*it)->getType()) {
707 if ((*it)->isEnabled()) {
708 x = (*it)->getResistance();
710 Node* nextNONremovableNode = (*it)->getTheOtherNode(node);
711 Element* nextSerialResistor = *it;
713 nextSerialResistor = nextNONremovableNode->
getAnOtherElement(nextSerialResistor);
715 nextNONremovableNode = nextSerialResistor->
getTheOtherNode(nextNONremovableNode);
719 eqn[node->
getId()] += x;
721 if (!nextNONremovableNode->
isGround()) {
722 eqn[nextNONremovableNode->
getId()] -= x;
727 if ((*it)->isEnabled()) {
729 if ((*it)->getPosNode() == node) {
732 x = (*it)->getPowerWanted() /
voltageSources->front()->getVoltage();
740 if ((*it)->getPosNode() == node) {
745 eqn[(*it)->getId()] += x;
747 removable_ids->push_back((*it)->getId());
766 for (std::vector<Node*>::iterator it =
nodes->begin(); it !=
nodes->end(); it++) {
768 if ((*it)->getElements()->size() == 2 && !(*it)->isGround()) {
770 (*it)->setRemovability(
true);
771 for (std::vector<Element*>::iterator it2 = (*it)->getElements()->begin(); it2 != (*it)->getElements()->end(); it2++) {
773 (*it)->setRemovability(
false);
777 if ((*it)->isRemovable()) {
779 removable_ids->push_back((*it)->getId());
782 (*it)->setRemovability(
false);
786 std::sort(removable_ids->begin(), removable_ids->end(), std::less<int>());
798 WRITE_WARNING(
TL(
"Trying to add resistor element into the overhead wire circuit with resistance < 1e-6. "))
800 WRITE_ERROR(
TL(
"Trying to add resistor element into the overhead wire circuit with resistance < 0. "))
809 std::cout <<
"The element: '" + name +
"' already exists.";
813 e =
new Element(name, et, value);
838 this->
elements->erase(std::remove(this->
elements->begin(), this->elements->end(), element), this->elements->end());
845 if (voltageSource->getNegNode() == unusedNode) {
846 voltageSource->setNegNode(newNode);
850 if (voltageSource->getPosNode() == unusedNode) {
851 voltageSource->setPosNode(newNode);
857 if (element->getNegNode() == unusedNode) {
858 element->setNegNode(newNode);
862 if (element->getPosNode() == unusedNode) {
863 element->setPosNode(newNode);
874 if (unusedNode->
getId() != modLastId) {
876 if (node_last !=
nullptr) {
880 if (elem_last !=
nullptr) {
883 WRITE_ERROR(
TL(
"The element or node with the last Id was not found in the circuit!"));
893 for (std::vector<Element*>::iterator it =
elements->begin(); it !=
elements->end(); it++) {
895 (*it)->setEnabled(
true);
900 (*it)->setEnabled(
true);
907 for (std::vector<Node*>::iterator it =
nodes->begin(); it !=
nodes->end(); it++) {
908 if ((*it)->getNumOfElements() < 2) {
910 if ((*it)->getNumOfElements() < 1) {
917 if ((*it)->getPosNode() ==
nullptr || (*it)->getNegNode() ==
nullptr) {
919 WRITE_ERRORF(
TL(
"Circuit Voltage Source '%' is connected to less than two nodes, please adjust the definition of the section (with substation '%')."), (*it)->getName(), substationId);
924 for (std::vector<Element*>::iterator it =
elements->begin(); it !=
elements->end(); it++) {
925 if ((*it)->getPosNode() ==
nullptr || (*it)->getNegNode() ==
nullptr) {
927 WRITE_ERRORF(
TL(
"Circuit Element '%' is connected to less than two nodes, please adjust the definition of the section (with substation '%')."), (*it)->getName(), substationId);
934 bool* nodesVisited =
new bool[num];
935 for (
int i = 0; i < num; i++) {
936 nodesVisited[i] =
false;
940 if (!
getNode(-1)->isGround()) {
942 WRITE_ERRORF(
TL(
"Circuit Node with id '-1' is not the grounded, please adjust the definition of the section (with substation '%')."), substationId);
944 std::vector<Node*>* queue =
new std::vector<Node*>(0);
945 Node* node =
nullptr;
946 Node* neigboringNode =
nullptr;
950 queue->push_back(node);
952 while (!queue->empty()) {
953 node = queue->back();
955 if (!nodesVisited[node->
getId()]) {
956 nodesVisited[node->
getId()] =
true;
958 neigboringNode = (*it)->getTheOtherNode(node);
960 queue->push_back(neigboringNode);
963 nodesVisited[(*it)->getId()] = 1;
966 WRITE_ERRORF(
TL(
"A Circuit Resistor Element connects the ground, please adjust the definition of the section (with substation '%')."), substationId);
972 for (
int i = 0; i < num; i++) {
973 if (nodesVisited[i] == 0) {
975 WRITE_WARNINGF(
TL(
"Circuit Node or Voltage Source with internal id '%' has been not visited during checking of the circuit. The circuit is disconnected, please adjust the definition of the section (with substation '%')."),
toString(i), substationId);
static std::mutex circuit_lock
#define WRITE_WARNINGF(...)
#define WRITE_ERRORF(...)
#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)
@ CURRENT_SOURCE_traction_wire
@ VOLTAGE_SOURCE_traction_wire
static bool gOverheadWireCurrentLimits
void setVoltage(double volt)
void setGround(bool isground)
void setNumMatrixCol(int num)
Element * getAnOtherElement(Element *element)
void addElement(Element *element)
void eraseElement(Element *element)
std::vector< Element * > * getElements()