diff --git a/README.md b/README.md index 5b869f7..033ef60 100644 --- a/README.md +++ b/README.md @@ -49,7 +49,6 @@ Input node data file and edge data file should be of CSV format. Definition of e * Load power (active power, PQ and PV nodes) * Load power (reactive power, PQ nodes) * (When three-phase short circuit calculation is enabled) Generator admittance -* (When three-phase short circuit calculation is enabled) Generator sub-transient voltage * Node type (0 - swing node, 1- PQ node, 2 - PV node) #### 2.2.2 Edge data file diff --git a/examples/3-gen-9-nodes/nodes.csv b/examples/3-gen-9-nodes/nodes.csv index 9692c50..26cb951 100644 --- a/examples/3-gen-9-nodes/nodes.csv +++ b/examples/3-gen-9-nodes/nodes.csv @@ -1,10 +1,10 @@ -U,Generator,P,Q,Xd,E,type -1.04,0,0,0,0.3,1.137,0 -1.025,1.63,0,0,0.3,1.211,2 -1.025,0.85,0,0,0.3,1.043,2 -1,0,0,0,0,0,1 -1,0,1.25,0.5,0,0,1 -1,0,0.9,0.3,0,0,1 -1,0,0,0,0,0,1 -1,0,1,0.35,0,0,1 -1,0,0,0,0,0,1 \ No newline at end of file +U,Generator,P,Q,Xd,type +1.04,0,0,0,0.3,0 +1.025,1.63,0,0,0.3,2 +1.025,0.85,0,0,0.3,2 +1,0,0,0,0,1 +1,0,1.25,0.5,0,1 +1,0,0.9,0.3,0,1 +1,0,0,0,0,1 +1,0,1,0.35,0,1 +1,0,0,0,0,1 \ No newline at end of file diff --git a/src/calc.cpp b/src/calc.cpp index f471b8a..c203b12 100644 --- a/src/calc.cpp +++ b/src/calc.cpp @@ -96,12 +96,12 @@ namespace flow bool verbose, double epsilon, bool short_circuit, bool ignore_load, unsigned short_circuit_node, const std::complex& z_f) { - if (nodes.n_cols != (short_circuit ? 7 : 5) || edges.n_cols != 6) + if (nodes.n_cols != (short_circuit ? 6 : 5) || edges.n_cols != 6) writer::error("Bad input matrix format."); short_circuit_ = short_circuit; nodes.each_row([this](const arma::rowvec& row) { - const auto type_val = static_cast(row[short_circuit_ ? 6 : 4]); + const auto type_val = static_cast(row[short_circuit_ ? 5 : 4]); auto type = node_data::swing; if (type_val == 1) { type = node_data::pq; @@ -113,11 +113,11 @@ namespace flow writer::error("Bad node type."); if (short_circuit_) { nodes_.push_back({ - num_nodes_++, row[0], row[1], row[2], row[3], row[4], row[5], type + num_nodes_++, row[0], row[1], row[2], row[3], row[4], type }); } else { nodes_.push_back({ - num_nodes_++, row[0], row[1], row[2], row[3], 0, 0, type + num_nodes_++, row[0], row[1], row[2], row[3], 0, type }); } }); @@ -143,38 +143,40 @@ namespace flow verbose_ = verbose; epsilon_ = epsilon; ignore_load_ = ignore_load; - short_circuit_node_ = short_circuit_node - 1; + if (short_circuit_node < 1 || short_circuit_node > num_nodes_) + writer::error("Bad node ID for short circuit calculation."); + short_circuit_node_ = node_offset(short_circuit_node - 1); z_f_ = z_f; } std::pair calc::node_admittance() { n_adm_.zeros(num_nodes_, num_nodes_); - arma::cx_mat adm_orig(num_nodes_, num_nodes_, arma::fill::zeros); + n_adm_orig_.zeros(num_nodes_, num_nodes_); for (auto&& edge : edges_) { - const auto impedance = edge.impedance(); + const auto admittance = edge.admittance(); const auto m = node_offset(edge.m); const auto n = node_offset(edge.n); // Whether this edge has transformer. if (edge.k) { - n_adm_.at(m, m) = adm_orig.at(edge.m, edge.m) += impedance; - n_adm_.at(n, n) = adm_orig.at(edge.n, edge.n) += impedance / std::pow(edge.k, 2); - n_adm_.at(m, n) = adm_orig.at(edge.m, edge.n) -= impedance / edge.k; - n_adm_.at(n, m) = adm_orig.at(edge.n, edge.m) -= impedance / edge.k; + n_adm_.at(m, m) = n_adm_orig_.at(edge.m, edge.m) += admittance; + n_adm_.at(n, n) = n_adm_orig_.at(edge.n, edge.n) += admittance / std::pow(edge.k, 2); + n_adm_.at(m, n) = n_adm_orig_.at(edge.m, edge.n) -= admittance / edge.k; + n_adm_.at(n, m) = n_adm_orig_.at(edge.n, edge.m) -= admittance / edge.k; } else { - const auto delta_diag = impedance + edge.grounding_admittance(); - n_adm_.at(m, m) = adm_orig.at(edge.m, edge.m) += delta_diag; - n_adm_.at(n, n) = adm_orig.at(edge.n, edge.n) += delta_diag; - n_adm_.at(m, n) = adm_orig.at(edge.m, edge.n) -= impedance; - n_adm_.at(n, m) = adm_orig.at(edge.n, edge.m) -= impedance; + const auto delta_diag = admittance + edge.grounding_admittance(); + n_adm_.at(m, m) = n_adm_orig_.at(edge.m, edge.m) += delta_diag; + n_adm_.at(n, n) = n_adm_orig_.at(edge.n, edge.n) += delta_diag; + n_adm_.at(m, n) = n_adm_orig_.at(edge.m, edge.n) -= admittance; + n_adm_.at(n, m) = n_adm_orig_.at(edge.n, edge.m) -= admittance; } adj_.at(m, n) = 1; adj_.at(n, m) = 1; } n_adm_g_ = arma::real(n_adm_); n_adm_b_ = arma::imag(n_adm_); - const auto n_adm_orig_g = arma::real(adm_orig); - const auto n_adm_orig_b = arma::imag(adm_orig); + const auto n_adm_orig_g = arma::real(n_adm_orig_); + const auto n_adm_orig_b = arma::imag(n_adm_orig_); if (verbose_) { writer::println("Real part of node admittance matrix:"); writer::print_mat(n_adm_orig_g); @@ -187,15 +189,17 @@ namespace flow std::pair calc::node_impedance() { n_imp_.zeros(num_nodes_, num_nodes_); - auto adm = n_adm_, adm_orig = n_adm_; + auto adm = n_adm_, adm_orig = n_adm_orig_; auto i = 0U; for (auto&& node : nodes_) { - const auto j = node_offset(i); if (node.type == node_data::pq) { - const auto impedance = std::complex(p_[i], -q_[i]) / std::pow(v_[i], 2); - adm.at(j, j) = adm_orig.at(i, i) += impedance; + if (!ignore_load_) { + // Note that we should use P(LD) and Q(LD). + const auto impedance = std::complex(-init_p_[i], init_q_[i]) / std::pow(v_[i], 2); + adm.at(i, i) = adm_orig.at(node.id, node.id) += impedance; + } } else { - adm.at(j, j) = adm_orig.at(i, i) -= std::complex(0, 1 / node.x_d); + adm.at(i, i) = adm_orig.at(node.id, node.id) -= std::complex(0, 1 / node.x_d); } ++i; } @@ -389,13 +393,11 @@ namespace flow std::complex calc::short_circuit_current() { const auto n = short_circuit_node_; - if (n > num_nodes_) - writer::error("Bad node ID for short circuit calculation."); if (ignore_load_) i_f_ = 1; else i_f_ = { e_[n], f_[n] }; - i_f_ /= std::complex(n_imp_b_.at(n, n), n_imp_g_.at(n, n)) + z_f_; + i_f_ /= std::complex(n_imp_g_.at(n, n), n_imp_b_.at(n, n)) + z_f_; return i_f_; } @@ -406,13 +408,17 @@ namespace flow const auto n = short_circuit_node_; vec_elem_foreach(u_f_, [n, this](auto& elem, auto row) { - std::complex u_f(1); - if (!ignore_load_) + std::complex u_f(1), u_i(1); + if (!ignore_load_) { u_f = { e_[row], f_[row] }; - elem = u_f - cx(n_imp_b_.at(row, n), n_imp_g_.at(row, n)) * - (cx(e_[n], f_[n]) / (cx(n_imp_b_.at(n, n), n_imp_g_.at(n, n)) + z_f_)); - if (approx_zero(elem.real()) && approx_zero(elem.imag())) - elem = { }; + u_i = { e_[n], f_[n] }; + } + elem = u_f - cx(n_imp_g_.at(row, n), n_imp_b_.at(row, n)) * + (u_i / (cx(n_imp_g_.at(n, n), n_imp_b_.at(n, n)) + z_f_)); + if (approx_zero(elem.real())) + elem.real(0); + if (approx_zero(elem.imag())) + elem.imag(0); }); arma::cx_colvec u_f_orig(num_nodes_); vec_elem_foreach(u_f_, [&u_f_orig, this](auto&& elem, auto row) @@ -421,10 +427,10 @@ namespace flow }); if (verbose_) { writer::println("Short circuit node voltage:"); - for (auto&& elem : u_f_) + for (auto&& elem : u_f_orig) writer::print_complex("", elem); } - return join_rows(arma::real(u_f_orig), arma::imag(u_f_orig)); + return join_rows(arma::real(u_f_), arma::imag(u_f_)); } arma::mat calc::short_circuit_edge_current() @@ -435,16 +441,16 @@ namespace flow writer::println("Short circuit edge current:"); for (auto&& edge : edges_) { - const auto impedance = edge.impedance(); + const auto admittance = edge.admittance(); const auto m = node_offset(edge.m); const auto n = node_offset(edge.n); - std::complex z; + std::complex y; if (edge.k) { - z = impedance * (edge.k * edge.k - edge.k + 1) / (edge.k * edge.k); + y = admittance * (edge.k * edge.k - edge.k + 1) / (edge.k * edge.k); } else { - z = impedance + edge.grounding_admittance() * 2.0; + y = admittance + edge.grounding_admittance() * 2.0; } - edge_current[i] = (u_f_[m] - u_f_[n] / (edge.k ? edge.k : 1)) / z; + edge_current[i] = (u_f_[m] - u_f_[n] / (edge.k ? edge.k : 1)) * y; if (verbose_) writer::print_complex(std::to_string(edge.m + 1) + ',' + std::to_string(edge.n + 1) + ": ", edge_current[i]); diff --git a/src/calc.hpp b/src/calc.hpp index aa3e1ab..d9475dd 100644 --- a/src/calc.hpp +++ b/src/calc.hpp @@ -35,9 +35,6 @@ namespace flow /// Generator admittance double x_d; - /// Sub-transient voltage of generator - double e; - /// Node type. enum node_type { pq, pv, swing @@ -63,16 +60,26 @@ namespace flow double k; /** - * Get impedance of edge. + * Get admittance of edge. * - * @return Impedance (complex). + * @return Admittance (complex). */ - std::complex impedance() const + std::complex admittance() const { const auto deno = r * r + x * x; return { r / deno, -x / deno }; } + /** + * Get impedance of edge. + * + * @return Admittance (complex). + */ + std::complex impedance() const + { + return { r, x }; + } + /** * Get grounding admittance. * @@ -102,6 +109,9 @@ namespace flow /// Node admittance matrix. arma::cx_mat n_adm_; + /// Node admittance matrix in original node order. + arma::cx_mat n_adm_orig_; + /// Node admittance matrix. arma::mat n_adm_g_, n_adm_b_;