30 #include "Eigen/Dense"
31 #include "Eigen/Sparse"
32 #include "Eigen/Geometry"
42 if (getNode(name) !=
nullptr) {
44 std::cout <<
"The node '" + name +
"' already exists." << std::endl;
48 if (nodes->size() == 0) {
51 Node* tNode =
new Node(name, this->lastId);
57 this->nodes->push_back(tNode);
64 this->nodes->erase(std::remove(this->nodes->begin(), this->nodes->end(), node), this->nodes->end());
69 Element* tElement = getElement(name);
70 if (tElement ==
nullptr) {
77 Element* tElement = getElement(name);
78 if (tElement ==
nullptr) {
79 Node* node = getNode(name);
80 if (node !=
nullptr) {
91 Element* tElement = getElement(name);
92 if (tElement ==
nullptr) {
100 for (
auto&& it : *this->nodes) {
101 if (it->getName() == name) {
109 for (vector<Node*>::iterator it = this->nodes->begin(); it != nodes->end(); it++) {
110 if ((*it)->getId() ==
id) {
118 for (vector<Element*>::iterator it = this->elements->begin(); it != elements->end(); it++) {
119 if ((*it)->getName() == name) {
123 for (vector<Element*>::iterator it = this->voltageSources->begin(); it != voltageSources->end(); it++) {
124 if ((*it)->getName() == name) {
132 for (vector<Element*>::iterator it = this->elements->begin(); it != elements->end(); it++) {
133 if ((*it)->getId() ==
id) {
137 for (vector<Element*>::iterator it = this->voltageSources->begin(); it != voltageSources->end(); it++) {
138 if ((*it)->getId() ==
id) {
146 for (vector<Element*>::iterator it = this->voltageSources->begin(); it != voltageSources->end(); it++) {
147 if ((*it)->getId() ==
id) {
155 vector<Element*>* vsources =
new vector<Element*>(0);
156 for (vector<Element*>::iterator it = this->elements->begin(); it != elements->end(); it++) {
157 if ((*it)->getType() == Element::ElementType::CURRENT_SOURCE_traction_wire) {
159 vsources->push_back(*it);
174 void Circuit::removeColumn(Eigen::MatrixXd& matrix,
int colToRemove) {
175 const int numRows = (int)matrix.rows();
176 const int numCols = (int)matrix.cols() - 1;
178 if (colToRemove < numCols) {
179 matrix.block(0, colToRemove, numRows, numCols - colToRemove) = matrix.rightCols(numCols - colToRemove);
182 matrix.conservativeResize(numRows, numCols);
185 bool Circuit::solveEquationsNRmethod(
double* eqn,
double* vals, std::vector<int>* removable_ids) {
187 int numofcolumn = (int)voltageSources->size() + (int)nodes->size() - 1;
188 int numofeqs = numofcolumn - (int)removable_ids->size();
190 Eigen::MatrixXd A = Eigen::Map < Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> >(eqn, numofeqs, numofcolumn);
194 for (std::vector<int>::reverse_iterator it = removable_ids->rbegin(); it != removable_ids->rend(); ++it) {
195 id = (*it >= 0 ? *it : -(*it));
203 Node* tNode =
nullptr;
204 for (
int i = 0; i < numofcolumn; i++) {
206 if (tNode !=
nullptr)
212 if (j > numofcolumn - (
int) removable_ids->size()) {
213 WRITE_ERROR(
"Number of column deployment during circuit evaluation was unsuccessfull.");
220 tElem = getElement(i);
221 if (tElem !=
nullptr) {
223 if (j > numofcolumn - (
int) removable_ids->size()) {
224 WRITE_ERROR(
"Number of column deployment deployment during circuit evaluation was unsuccessfull.");
230 WRITE_ERROR(
"Number of column deployment during circuit evaluation was unsuccessfull.");
233 Eigen::Map<Eigen::VectorXd> b(vals, numofeqs);
234 Eigen::VectorXd x = A.colPivHouseholderQr().solve(b);
238 Eigen::MatrixXd J = A;
240 int max_iter_of_NR = 10;
245 std::vector<double> alpha_notSolution;
246 double alpha_res = 1e-2;
247 double* x_best =
new double[numofeqs];
249 for (
int i = 0; i < numofeqs; i++) {
252 if (x.maxCoeff() > 10e6 || x.minCoeff() < -10e6) {
257 for (
int i = 0; i < numofeqs; i++) {
269 for (
int i = 0; i < numofeqs - (int) voltageSources->size(); i++) {
275 for (
auto& node : *nodes) {
276 if (node->isGround() || node->isRemovable() || node->getNumMatrixRow() == -2) {
279 if (node->getNumMatrixRow() != i) {
280 WRITE_ERROR(
"wrongly assigned row of matrix A during solving the circuit");
283 for (
auto it_element = node->getElements()->begin(); it_element != node->getElements()->end(); it_element++) {
284 if ((*it_element)->getType() == Element::ElementType::CURRENT_SOURCE_traction_wire) {
285 if ((*it_element)->isEnabled()) {
287 if ((*it_element)->getPosNode()->getNumMatrixCol() == -1) {
288 diff_voltage = -x[(*it_element)->getNegNode()->getNumMatrixCol()];
289 }
else if ((*it_element)->getNegNode()->getNumMatrixCol() == -1) {
290 diff_voltage = x[(*it_element)->getPosNode()->getNumMatrixCol()];
292 diff_voltage = (x[(*it_element)->getPosNode()->getNumMatrixCol()] - x[(*it_element)->getNegNode()->getNumMatrixCol()]);
295 if ((*it_element)->getPosNode() == node) {
296 vals[i] -= alpha * (*it_element)->getPowerWanted() / diff_voltage;
297 (*it_element)->setCurrent(-alpha * (*it_element)->getPowerWanted() / diff_voltage);
298 if ((*it_element)->getPosNode()->getNumMatrixCol() != -1) {
299 J(i, (*it_element)->getPosNode()->getNumMatrixCol()) -= alpha * (*it_element)->getPowerWanted() / diff_voltage / diff_voltage;
301 if ((*it_element)->getNegNode()->getNumMatrixCol() != -1) {
302 J(i, (*it_element)->getNegNode()->getNumMatrixCol()) += alpha * (*it_element)->getPowerWanted() / diff_voltage / diff_voltage;
305 vals[i] += alpha * (*it_element)->getPowerWanted() / diff_voltage;
308 if ((*it_element)->getPosNode()->getNumMatrixCol() != -1) {
309 J(i, (*it_element)->getPosNode()->getNumMatrixCol()) += alpha * (*it_element)->getPowerWanted() / diff_voltage / diff_voltage;
311 if ((*it_element)->getNegNode()->getNumMatrixCol() != -1) {
312 J(i, (*it_element)->getNegNode()->getNumMatrixCol()) -= alpha * (*it_element)->getPowerWanted() / diff_voltage / diff_voltage;
322 Eigen::Map<Eigen::VectorXd> bb(vals, numofeqs);
324 if ((A * x - bb).norm() < 1e-6) {
326 for (
int ii = 0; ii < numofeqs; ii++) {
330 }
else if (iterNR == max_iter_of_NR) {
331 alpha_notSolution.push_back(alpha);
332 for (
int ii = 0; ii < numofeqs; ii++) {
338 dx = -J.colPivHouseholderQr().solve(A * x - bb);
343 if (alpha_notSolution.empty()) {
347 if ((alpha_notSolution.back() - alphaBest) < alpha_res) {
348 max_iter_of_NR = 2 * max_iter_of_NR;
349 alpha_res = alpha_res / 10;
350 if (alpha_res < 5e-5) {
353 alpha = alpha_notSolution.back();
354 alpha_notSolution.pop_back();
358 alpha = alphaBest + 0.5 * (alpha_notSolution.back() - alphaBest);
361 for (
int i = 0; i < numofeqs; i++) {
366 for (
auto& node : *nodes) {
367 if (node->isGround() || node->isRemovable() || node->getNumMatrixRow() == -2) {
370 if (node->getNumMatrixRow() != i) {
371 WRITE_ERROR(
"wrongly assigned row of matrix A during solving the circuit");
373 for (
auto it_element = node->getElements()->begin(); it_element != node->getElements()->end(); it_element++) {
374 if ((*it_element)->getType() == Element::ElementType::CURRENT_SOURCE_traction_wire) {
375 if ((*it_element)->isEnabled()) {
377 if ((*it_element)->getPosNode()->getNumMatrixCol() == -1) {
378 diff_voltage = -x_best[(*it_element)->getNegNode()->getNumMatrixCol()];
379 }
else if ((*it_element)->getNegNode()->getNumMatrixCol() == -1) {
380 diff_voltage = x_best[(*it_element)->getPosNode()->getNumMatrixCol()];
382 diff_voltage = (x_best[(*it_element)->getPosNode()->getNumMatrixCol()] - x_best[(*it_element)->getNegNode()->getNumMatrixCol()]);
385 if ((*it_element)->getPosNode() == node) {
386 (*it_element)->setCurrent(-alphaBest * (*it_element)->getPowerWanted() / diff_voltage);
402 int n = (int)(voltageSources->size() + nodes->size() - 1);
405 Node* tNode =
nullptr;
406 for (
int i = 0; i < n; i++) {
408 if (tNode !=
nullptr)
413 if (j > n - (
int) removable_ids->size()) {
414 WRITE_ERROR(
"Results deployment during circuit evaluation was unsuccessfull.");
421 tElem = getElement(i);
422 if (tElem !=
nullptr) {
424 if (j > n - (
int) removable_ids->size()) {
425 WRITE_ERROR(
"Results deployment during circuit evaluation was unsuccessfull.");
435 WRITE_ERROR(
"Results deployment during circuit evaluation was unsuccessfull.");
440 Node* nextNONremovableNode1 =
nullptr;
441 Node* nextNONremovableNode2 =
nullptr;
443 for (vector<Node*>::iterator it = nodes->begin(); it != nodes->end(); it++) {
444 if (!(*it)->isRemovable()) {
447 if ((*it)->getElements()->size() != 2) {
451 el1 = (*it)->getElements()->front();
452 el2 = (*it)->getElements()->back();
471 y = ((1 - x) * nextNONremovableNode1->
getVoltage()) + (x * nextNONremovableNode2->
getVoltage());
473 (*it)->setRemovability(
false);
493 nodes =
new vector<Node*>(0);
494 elements =
new vector<Element*>(0);
495 voltageSources =
new vector<Element*>(0);
501 bool Circuit::_solveNRmethod() {
502 double* eqn =
nullptr;
503 double* vals =
nullptr;
504 std::vector<int> removable_ids;
506 detectRemovableNodes(&removable_ids);
507 createEquationsNRmethod(eqn, vals, &removable_ids);
508 if (!solveEquationsNRmethod(eqn, vals, &removable_ids)) {
511 deployResults(vals, &removable_ids);
516 bool Circuit::solve() {
520 return this->_solveNRmethod();
523 bool Circuit::createEquationsNRmethod(
double*& eqs,
double*& vals, std::vector<int>* removable_ids) {
525 int n = (int)(voltageSources->size() + nodes->size() - 1);
526 int m = n - (int)(removable_ids->size() - voltageSources->size());
528 eqs =
new double[m * n];
529 vals =
new double[m];
531 for (
int i = 0; i < m; i++) {
533 for (
int j = 0; j < n; j++) {
539 for (vector<Node*>::iterator it = nodes->begin(); it != nodes->end(); it++) {
540 if ((*it)->isGround() || (*it)->isRemovable()) {
541 (*it)->setNumMatrixRow(-1);
544 bool noVoltageSource = createEquationNRmethod((*it), (eqs + n * i), vals[i], removable_ids);
546 if (noVoltageSource) {
547 (*it)->setNumMatrixRow(i);
550 (*it)->setNumMatrixRow(-2);
552 for (
int j = 0; j < n; j++) {
558 std::sort(removable_ids->begin(), removable_ids->end(), std::less<int>());
560 for (vector<Element*>::iterator it = voltageSources->begin(); it != voltageSources->end(); it++) {
561 createEquation((*it), (eqs + n * i), vals[i]);
568 bool Circuit::createEquation(
Element* vsource,
double* eqn,
double& val) {
583 bool Circuit::createEquationNRmethod(
Node* node,
double* eqn,
double& val, std::vector<int>* removable_ids) {
584 for (vector<Element*>::iterator it = node->
getElements()->begin(); it != node->
getElements()->end(); it++) {
586 switch ((*it)->getType()) {
587 case Element::ElementType::RESISTOR_traction_wire:
588 if ((*it)->isEnabled()) {
589 x = (*it)->getResistance();
590 Node* nextNONremovableNode = (*it)->getTheOtherNode(node);
591 Element* nextSerialResistor = *it;
593 nextSerialResistor = nextNONremovableNode->
getAnOtherElement(nextSerialResistor);
595 nextNONremovableNode = nextSerialResistor->
getTheOtherNode(nextNONremovableNode);
598 eqn[node->
getId()] += x;
599 if (!nextNONremovableNode->
isGround()) {
600 eqn[nextNONremovableNode->
getId()] -= x;
606 case Element::ElementType::CURRENT_SOURCE_traction_wire:
607 if ((*it)->isEnabled()) {
608 if ((*it)->getPosNode() == node) {
609 x = (*it)->getCurrent();
611 x = -(*it)->getCurrent();
618 case Element::ElementType::VOLTAGE_SOURCE_traction_wire:
619 if ((*it)->getPosNode() == node) {
624 eqn[(*it)->getId()] += x;
626 removable_ids->push_back((*it)->getId());
629 case Element::ElementType::ERROR_traction_wire:
639 for (vector<Node*>::iterator it = nodes->begin(); it != nodes->end(); it++) {
640 if ((*it)->getElements()->size() == 2 && !(*it)->isGround()) {
641 (*it)->setRemovability(
true);
642 for (vector<Element*>::iterator it2 = (*it)->getElements()->begin(); it2 != (*it)->getElements()->end(); it2++) {
643 if ((*it2)->getType() != Element::ElementType::RESISTOR_traction_wire) {
644 (*it)->setRemovability(
false);
648 if ((*it)->isRemovable()) {
649 removable_ids->push_back((*it)->getId());
652 (*it)->setRemovability(
false);
655 std::sort(removable_ids->begin(), removable_ids->end(), std::less<int>());
660 if ((et == Element::ElementType::RESISTOR_traction_wire && value <= 0) || et == Element::ElementType::ERROR_traction_wire) {
668 std::cout <<
"The element: '" + name +
"' already exists.";
672 e =
new Element(name, et, value);
673 if (e->
getType() == Element::ElementType::VOLTAGE_SOURCE_traction_wire) {
677 this->voltageSources->push_back(e);
681 this->elements->push_back(e);
697 this->elements->erase(std::remove(this->elements->begin(), this->elements->end(), element), this->elements->end());
703 for (
auto& voltageSource : *voltageSources) {
704 if (voltageSource->getNegNode() == unusedNode) {
705 voltageSource->setNegNode(newNode);
709 if (voltageSource->getPosNode() == unusedNode) {
710 voltageSource->setPosNode(newNode);
715 for (
auto& element : *elements) {
716 if (element->getNegNode() == unusedNode) {
717 element->setNegNode(newNode);
721 if (element->getPosNode() == unusedNode) {
722 element->setPosNode(newNode);
729 this->eraseNode(unusedNode);
732 int modLastId = this->getLastId() - 1;
733 if (unusedNode->
getId() != modLastId) {
734 Node* node_last = this->getNode(modLastId);
735 if (node_last !=
nullptr) {
738 Element* elem_last = this->getVoltageSource(modLastId);
739 if (elem_last !=
nullptr) {
742 WRITE_ERROR(
"The element or node with the last Id was not found in the circuit!");
747 this->descreaseLastId();
752 for (vector<Element*>::iterator it = elements->begin(); it != elements->end(); it++) {
753 if ((*it)->getType() != Element::ElementType::RESISTOR_traction_wire) {
754 (*it)->setEnabled(
true);
758 for (vector<Element*>::iterator it = voltageSources->begin(); it != voltageSources->end(); it++) {
759 (*it)->setEnabled(
true);
761 this->iscleaned =
true;
766 for (vector<Node*>::iterator it = nodes->begin(); it != nodes->end(); it++) {
767 if ((*it)->getNumOfElements() < 2) {
769 if ((*it)->getNumOfElements() < 1) {
775 for (vector<Element*>::iterator it = voltageSources->begin(); it != voltageSources->end(); it++) {
776 if ((*it)->getPosNode() ==
nullptr || (*it)->getNegNode() ==
nullptr) {
778 WRITE_ERROR(
"Circuit Voltage Source '" + (*it)->getName() +
"' is connected to less than two nodes, please adjust the definition of the section (with substation '" + substationId +
"').");
783 for (vector<Element*>::iterator it = elements->begin(); it != elements->end(); it++) {
784 if ((*it)->getPosNode() ==
nullptr || (*it)->getNegNode() ==
nullptr) {
786 WRITE_ERROR(
"Circuit Element '" + (*it)->getName() +
"' is connected to less than two nodes, please adjust the definition of the section (with substation '" + substationId +
"').");
792 int num = (int)nodes->size() + getNumVoltageSources() - 1;
793 bool* nodesVisited =
new bool[num];
794 for (
int i = 0; i < num; i++) {
795 nodesVisited[i] =
false;
799 if (!getNode(-1)->isGround()) {
801 WRITE_ERROR(
"Circuit Node with id '-1' is not the grounded, please adjust the definition of the section (with substation '" + substationId +
"').");
803 vector<Node*>* queue =
new vector<Node*>(0);
804 Node* node =
nullptr;
805 Node* neigboringNode =
nullptr;
807 nodesVisited[voltageSources->front()->getId()] = 1;
808 node = voltageSources->front()->getPosNode();
809 queue->push_back(node);
811 while (!queue->empty()) {
812 node = queue->back();
814 if (!nodesVisited[node->
getId()]) {
815 nodesVisited[node->
getId()] =
true;
817 neigboringNode = (*it)->getTheOtherNode(node);
819 queue->push_back(neigboringNode);
820 }
else if ((*it)->getType() == Element::ElementType::VOLTAGE_SOURCE_traction_wire) {
822 nodesVisited[(*it)->getId()] = 1;
823 }
else if ((*it)->getType() == Element::ElementType::RESISTOR_traction_wire) {
825 WRITE_ERROR(
"A Circuit Resistor Element connects the ground, please adjust the definition of the section (with substation '" + substationId +
"').");
831 for (
int i = 0; i < num; i++) {
832 if (nodesVisited[i] == 0) {
834 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 +
"').");
842 return (
int) voltageSources->size();
#define WRITE_WARNING(msg)
std::string toString(const T &t, std::streamsize accuracy=gPrecision)
double getCurrent(string name)
int getNumVoltageSources()
vector< Element * > * getCurrentSources()
Element * addElement(string name, double value, Node *pNode, Node *nNode, Element::ElementType et)
void eraseNode(Node *node)
void detectRemovableNodes(std::vector< int > *removable_ids)
void replaceAndDeleteNode(Node *unusedNode, Node *newNode)
Node * getNode(string name)
double getVoltage(string name)
Element * getVoltageSource(int id)
void deployResults(double *vals, std::vector< int > *removable_ids)
Element * getElement(string name)
bool checkCircuit(std::string substationId="")
Node * addNode(string name)
double getResistance(string name)
void eraseElement(Element *element)
void setNegNode(Node *node)
void setPosNode(Node *node)
Node * getTheOtherNode(Node *node)
vector< Element * > * getElements()
void setVoltage(double voltage)
void setGround(bool isground)
void setNumMatrixCol(int num)
Element * getAnOtherElement(Element *element)
void addElement(Element *element)
void eraseElement(Element *element)