archived
/
arma-flow
Archived
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Add three-phase short circuit calculation.

This commit is contained in:
CismonX 2018-05-31 11:16:44 +08:00
parent b922b3da13
commit 74913850a3
10 changed files with 328 additions and 63 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
README.md
View File

@ -1,13 +1,16 @@
# arma-flow
A simple power flow calculator using Newton's method, written in C++ and armadillo.
A simple power flow calculator using Newton's method, written in C++ using armadillo.
Additional support for three-phase short-circuit calculation.
## 1. Requirements
* A newer version of [Armadillo](http://arma.sourceforge.net/)
* Tested on 8.400.0
* The [Options](https://mulholland.xyz/docs/options/) library
* Compiler with supports C++17
* This library was recently renamed to Janus. In order to build correctly, you may need to modify some code.
* Compiler with C++17 support
## 2. Documentation
@ -21,17 +24,18 @@ A simple power flow calculator using Newton's method, written in C++ and armadil
* `-r` : Remove first line of input CSV files before parsing.
* `-i <max_iterations>` : Max number of iterations to be performed before aborting.
* `-a <accuracy>` : Max deviation to be tolerated.
* `-s <node_id>` : Calculate three-phase short circuit on specified node.
* `--ignore-load` : Ignore load current when calculating three-phase short circuit.
* `--tr <transition_impedance(real)>` : Transition impedance of three-phase short circuit(real part).
* `--ti <transition_impedance(imag)>` : Transition impedance of three-phase short circuit(imaginary part).
* `-v | --verbose` : Output more text to STDOUT.
For example:
```bash
# Building arma-flow
git clone https://github.com/CismonX/arma-flow.git
cd arma-flow
make -j4
# Testing arma-flow
./arma-flow -n examples/3-gen-9-nodes/nodes.csv -e examples/3-gen-9-nodes/edges.csv -r
./arma-flow -n examples/3-gen-9-nodes/nodes.csv \
-e examples/3-gen-9-nodes/edges.csv \
-r -s 4 --ignore-load --verbose
```
### 2.2 Input
@ -44,6 +48,8 @@ Input node data file and edge data file should be of CSV format. Definition of e
* Generator power (active power, PV nodes)
* 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
@ -67,9 +73,17 @@ The real part and imaginary part of node admittance matrix will be printed to "\
#### 2.3.2 Calculation result
The result of power flow calculation will be printed to "\<prefix\>result.csv". The definition of each column is given below:
The result of power flow calculation will be printed to "\<prefix\>flow.csv". The definition of each column is given below:
* Node voltage
* Phase angle of node voltage (Radian)
* Node power (active)
* Node power (reactive)
* Node power (reactive)
#### 2.3.3 Short circuit calculation
The node impedance matrix will be printed to "\<prefix\>node-impedance-real.csv" and "\<prefix\>node-impedance-imag.csv".
Short circuit current will be printed directly to STDOUT.
Node voltage after short circuit will be printed to "\<prefix\>short-circuit-voltage.csv", and edge current "\<prefix\>short-circuit-edge-current.csv".

20
examples/3-gen-9-nodes/nodes.csv
View File

@ -1,10 +1,10 @@
U,Generator,P,Q,type
1.04,0,0,0,0
1.025,1.63,0,0,2
1.025,0.85,0,0,2
1,0,0,0,1
1,0,1.25,0.5,1
1,0,0.9,0.3,1
1,0,0,0,1
1,0,1,0.35,1
1,0,0,0,1
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
1 U Generator P Q Xd E type
2 1.04 0 0 0 0.3 1.137 0
3 1.025 1.63 0 0 0.3 1.211 2
4 1.025 0.85 0 0 0.3 1.043 2
5 1 0 0 0 0 0 1
6 1 0 1.25 0.5 0 0 1
7 1 0 0.9 0.3 0 0 1
8 1 0 0 0 0 0 1
9 1 0 1 0.35 0 0 1
10 1 0 0 0 0 0 1

30
src/args.cpp
View File

@ -12,8 +12,9 @@ namespace flow
"A simple power flow calculator using Newton's method.\n"
"usage: arma-flow [--version] [-h | --help] [-o <output_file_prefix>]\n"
" -n <node_data_file> -e <edge_data_file> [-r]\n"
" [-i <max_iterations>] [-a <accuracy>]\n"
" [-v | --verbose]",
" [-i <max_iterations>] [-a <accuracy>] [-v | --verbose]\n"
" [-s <node_id>] [--ignore-load]\n"
" [--tr <transition_impedance(real)>] [--ti <transition_impedance(imag)>]",
"arma-flow version 0.0.1")
{
arg_parser_.newString("o", "result-");
@ -22,6 +23,10 @@ namespace flow
arg_parser_.newFlag("r");
arg_parser_.newInt("i", 100);
arg_parser_.newDouble("a", 0.00001);
arg_parser_.newInt("s");
arg_parser_.newFlag("ignore-load");
arg_parser_.newDouble("tr", 0);
arg_parser_.newDouble("ti", 0);
arg_parser_.newFlag("verbose v");
}
@ -62,6 +67,27 @@ namespace flow
return arg_parser_.found("a");
}
bool args::short_circuit(unsigned& node)
{
if (!arg_parser_.found("s"))
return false;
node = arg_parser_.getInt("s");
return true;
}
bool args::ignore_load()
{
return arg_parser_.getFlag("ignore-load");
}
bool args::transition_impedance(std::complex<double>& z_f)
{
if (!arg_parser_.found("tr") && !arg_parser_.found("ti"))
return false;
z_f = { arg_parser_.getDouble("tr"), arg_parser_.getDouble("ti") };
return true;
}
bool args::verbose()
{
return arg_parser_.getFlag("v");

20
src/args.hpp
View File

@ -6,6 +6,7 @@
#pragma once
#include <complex>
#include <options.h>
namespace flow
@ -49,8 +50,6 @@ namespace flow
/**
* Check whether to remove the first line from input files.
*
* @return Whether argument is provided.
*/
bool remove_first_line();
@ -68,6 +67,23 @@ namespace flow
*/
bool accuracy(double& epsilon);
/**
* Calculate three-phase short circuit on specified node.
*
* @return Whether argument is provided.
*/
bool short_circuit(unsigned& node);
/**
* Check whether to ignore load current when calculating short circuit.
*/
bool ignore_load();
/**
* Get transition impedance of three-phase short circuit.
*/
bool transition_impedance(std::complex<double>& z_f);
/**
* Check whether to enable verbose output.
*

153
src/calc.cpp
View File

@ -93,13 +93,15 @@ namespace flow
}
void calc::init(const arma::mat& nodes, const arma::mat& edges,
bool verbose, double epsilon)
bool verbose, double epsilon, bool short_circuit, bool ignore_load,
unsigned short_circuit_node, const std::complex<double>& z_f)
{
if (nodes.n_cols != 5 || edges.n_cols != 6)
if (nodes.n_cols != (short_circuit ? 7 : 5) || edges.n_cols != 6)
writer::error("Bad input matrix format.");
short_circuit_ = short_circuit;
nodes.each_row([this](const arma::rowvec& row)
{
auto type_val = static_cast<unsigned>(row[4]);
const auto type_val = static_cast<unsigned>(row[short_circuit_ ? 6 : 4]);
auto type = node_data::swing;
if (type_val == 1) {
type = node_data::pq;
@ -109,9 +111,15 @@ namespace flow
++num_pv_;
} else if (type_val != 0)
writer::error("Bad node type.");
nodes_.push_back({
num_nodes_++, row[0], row[1], row[2], row[3], type
});
if (short_circuit_) {
nodes_.push_back({
num_nodes_++, row[0], row[1], row[2], row[3], row[4], row[5], type
});
} else {
nodes_.push_back({
num_nodes_++, row[0], row[1], row[2], row[3], 0, 0, type
});
}
});
// Nodes should be sorted, PQ nodes should be followed by PV nodes,
// while swing node be the last.
@ -124,8 +132,8 @@ namespace flow
writer::error("Only one swing node should exist.");
edges.each_row([this](const arma::rowvec& row)
{
auto n1 = static_cast<unsigned>(row[0]) - 1;
auto n2 = static_cast<unsigned>(row[1]) - 1;
const auto n1 = static_cast<unsigned>(row[0]) - 1;
const auto n2 = static_cast<unsigned>(row[1]) - 1;
if (n1 >= num_nodes_ || n2 >= num_nodes_)
writer::error("Bad node offset.");
edges_.push_back({
@ -134,11 +142,14 @@ namespace flow
});
verbose_ = verbose;
epsilon_ = epsilon;
ignore_load_ = ignore_load;
short_circuit_node_ = short_circuit_node - 1;
z_f_ = z_f;
}
std::pair<arma::mat, arma::mat> calc::node_admittance()
{
arma::cx_mat node_adm_cplx(num_nodes_, num_nodes_, arma::fill::zeros);
n_adm_.zeros(num_nodes_, num_nodes_);
arma::cx_mat adm_orig(num_nodes_, num_nodes_, arma::fill::zeros);
for (auto&& edge : edges_) {
const auto impedance = edge.impedance();
@ -146,22 +157,22 @@ namespace flow
const auto n = node_offset(edge.n);
// Whether this edge has transformer.
if (edge.k) {
node_adm_cplx.at(m, m) += adm_orig(edge.m, edge.m) = impedance;
node_adm_cplx.at(n, n) += adm_orig(edge.n, edge.n) = impedance / std::pow(edge.k, 2);
node_adm_cplx.at(m, n) -= adm_orig(edge.m, edge.n) = impedance / edge.k;
node_adm_cplx.at(n, m) -= adm_orig(edge.n, edge.m) = impedance / 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;
} else {
const auto delta_diag = impedance + edge.grounding_admittance();
node_adm_cplx.at(m, m) += adm_orig(edge.m, edge.m) = delta_diag;
node_adm_cplx.at(n, n) += adm_orig(edge.n, edge.n) = delta_diag;
node_adm_cplx.at(m, n) -= adm_orig(edge.m, edge.n) = impedance;
node_adm_cplx.at(n, m) -= adm_orig(edge.n, edge.m) = impedance;
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;
}
adj_.at(m, n) = 1;
adj_.at(n, m) = 1;
}
n_adm_g_ = arma::real(node_adm_cplx);
n_adm_b_ = arma::imag(node_adm_cplx);
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);
if (verbose_) {
@ -173,6 +184,36 @@ namespace flow
return { n_adm_orig_g, n_adm_orig_b };
}
std::pair<arma::mat, arma::mat> calc::node_impedance()
{
n_imp_.zeros(num_nodes_, num_nodes_);
auto adm = n_adm_, adm_orig = n_adm_;
auto i = 0U;
for (auto&& node : nodes_) {
const auto j = node_offset(i);
if (node.type == node_data::pq) {
const auto impedance = std::complex<double>(p_[i], -q_[i]) / std::pow(v_[i], 2);
adm.at(j, j) = adm_orig.at(i, i) += impedance;
} else {
adm.at(j, j) = adm_orig.at(i, i) -= std::complex<double>(0, 1 / node.x_d);
}
++i;
}
n_imp_ = adm.i();
n_imp_g_ = arma::real(n_imp_);
n_imp_b_ = arma::imag(n_imp_);
const auto imp_orig = adm_orig.i();
const auto n_imp_orig_g = arma::real(imp_orig);
const auto n_imp_orig_b = arma::imag(imp_orig);
if (verbose_) {
writer::println("Real part of node impedance matrix:");
writer::print_mat(n_imp_orig_g);
writer::println("Imaginary part of node impedance matrix:");
writer::print_mat(n_imp_orig_b);
}
return { n_imp_orig_g, n_imp_orig_b };
}
void calc::iterate_init()
{
init_p_.zeros(num_nodes_ - 1);
@ -316,8 +357,8 @@ namespace flow
for (auto&& elem : q_)
if (approx_zero(elem))
elem = 0;
arma::colvec v(num_nodes_);
vec_elem_foreach(v, [this](auto& elem, auto col)
v_.zeros(num_nodes_);
vec_elem_foreach(v_, [this](auto& elem, auto col)
{
elem = std::sqrt(std::pow(e_[col], 2) + std::pow(f_[col], 2));
if (approx_zero(elem))
@ -335,7 +376,8 @@ namespace flow
{
elem = nodes_[row].id;
});
arma::mat retval = join_rows(node_id, join_rows(join_rows(v, theta), join_rows(p_, q_)));
arma::mat retval = join_rows(node_id, join_rows(join_rows(v_, theta), join_rows(p_, q_)));
// We shall preserve the original node sequence.
arma::mat sorted_retval(num_nodes_, 4);
retval.each_row([&sorted_retval](const arma::rowvec& row)
{
@ -343,4 +385,71 @@ namespace flow
});
return sorted_retval;
}
std::complex<double> 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<double>(n_imp_b_.at(n, n), n_imp_g_.at(n, n)) + z_f_;
return i_f_;
}
arma::mat calc::short_circuit_voltage()
{
u_f_.resize(num_nodes_);
using cx = std::complex<double>;
const auto n = short_circuit_node_;
vec_elem_foreach(u_f_, [n, this](auto& elem, auto row)
{
std::complex<double> u_f(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 = { };
});
arma::cx_colvec u_f_orig(num_nodes_);
vec_elem_foreach(u_f_, [&u_f_orig, this](auto&& elem, auto row)
{
u_f_orig[nodes_[row].id] = elem;
});
if (verbose_) {
writer::println("Short circuit node voltage:");
for (auto&& elem : u_f_)
writer::print_complex("", elem);
}
return join_rows(arma::real(u_f_orig), arma::imag(u_f_orig));
}
arma::mat calc::short_circuit_edge_current()
{
arma::cx_vec edge_current(edges_.size());
auto i = 0U;
if (verbose_)
writer::println("Short circuit edge current:");
for (auto&& edge : edges_)
{
const auto impedance = edge.impedance();
const auto m = node_offset(edge.m);
const auto n = node_offset(edge.n);
std::complex<double> z;
if (edge.k) {
z = impedance * (edge.k * edge.k - edge.k + 1) / (edge.k * edge.k);
} else {
z = impedance + edge.grounding_admittance() * 2.0;
}
edge_current[i] = (u_f_[m] - u_f_[n] / (edge.k ? edge.k : 1)) / z;
if (verbose_)
writer::print_complex(std::to_string(edge.m + 1) + ',' +
std::to_string(edge.n + 1) + ": ", edge_current[i]);
++i;
}
return join_rows(arma::real(edge_current), arma::imag(edge_current));
}
}

68
src/calc.hpp
View File

@ -32,6 +32,12 @@ namespace flow
/// Load (reactive power)
double q;
/// Generator admittance
double x_d;
/// Sub-transient voltage of generator
double e;
/// Node type.
enum node_type {
pq, pv, swing
@ -93,9 +99,18 @@ namespace flow
/// Adjacency matrix of nodes.
arma::uchar_mat adj_;
/// Node admittance matrix.
arma::cx_mat n_adm_;
/// Node admittance matrix.
arma::mat n_adm_g_, n_adm_b_;
/// Node impedance matrix.
arma::cx_mat n_imp_;
/// Node impedance matrix.
arma::mat n_imp_g_, n_imp_b_;
/// Given values of power and voltage.
arma::colvec init_p_, init_q_, init_v_;
@ -117,8 +132,33 @@ namespace flow
/// Power vector of nodes.
arma::colvec p_, q_;
/// Voltage of nodes.
arma::colvec v_;
/// Whether to calculate short circuit.
bool short_circuit_;
/// Whether to ignore load current when calculating short circuit.
bool ignore_load_;
/// Node ID of three-phase short circuit.
unsigned short_circuit_node_;
/// Transition impedance of node.
std::complex<double> z_f_;
/// Value of short circuit current.
std::complex<double> i_f_;
/// Vector of short circuit voltage.
arma::cx_colvec u_f_;
arma::mat i_real_;
arma::mat i_imag_;
/// Whether verbose output is enabled.
bool verbose_ = false;
bool verbose_;
/// Max deviation to be tolerated.
double epsilon_ = 0;
@ -228,13 +268,19 @@ namespace flow
* Initialize.
*/
void init(const arma::mat& nodes, const arma::mat& edges,
bool verbose, double epsilon);
bool verbose, double epsilon, bool short_circuit, bool ignore_load,
unsigned short_circuit_node, const std::complex<double>& z_f);
/**
* Calculate node admittance.
*/
std::pair<arma::mat, arma::mat> node_admittance();
/**
* Calculate node impedance.
*/
std::pair<arma::mat, arma::mat> node_impedance();
/**
* Initialize iteration.
*/
@ -253,9 +299,23 @@ namespace flow
double get_max() const;
/**
* Get result of calculation.
*
* Get result of power flow calculation.
*/
arma::mat result();
/**
* Get current of three-phase short circuit.
*/
std::complex<double> short_circuit_current();
/**
* Get node voltage of short circuit.
*/
arma::mat short_circuit_voltage();
/**
* Get edge current of short circuit.
*/
arma::mat short_circuit_edge_current();
};
}

37
src/executor.cpp
View File

@ -10,14 +10,15 @@
namespace flow
{
executor::executor() : factory_(factory::get()) {}
void executor::execute(int argc, char** argv) const
{
// Get components.
auto factory = factory::get();
auto args = factory->get_args();
auto input = factory->get_reader();
auto calc = factory->get_calc();
auto writer = factory->get_writer();
auto args = factory_->get_args();
auto input = factory_->get_reader();
auto calc = factory_->get_calc();
auto writer = factory_->get_writer();
// Parse options.
args->parse(argc, argv);
@ -48,9 +49,16 @@ namespace flow
if (!args->output_file_path(output_path) && verbose)
writer::notice("Output file path not specified. Defaulted to result-*.csv.");
writer->set_output_path_prefix(output_path);
unsigned short_circuit_node;
const auto short_circuit = args->short_circuit(short_circuit_node);
std::complex<double> transition_impedance;
if (short_circuit && !args->transition_impedance(transition_impedance) && verbose)
writer::notice("Transition impedance not specified, Defaulted to 0.");
const auto ignore_load = args->ignore_load();
// Initialize calculation.
calc->init(nodes, edges, verbose, epsilon);
calc->init(nodes, edges, verbose, epsilon, short_circuit, ignore_load,
short_circuit_node, transition_impedance);
const auto admittance = calc->node_admittance();
writer->to_csv_file("node-admittance-real.csv", admittance.first);
writer->to_csv_file("node-admittance-imag.csv", admittance.second);
@ -68,9 +76,24 @@ namespace flow
writer::print_mat(result);
}
writer->to_csv_file("flow.csv", result, "V,theta,P,Q");
return;
goto flow_done;
}
} while (num_iterations < max);
writer::error("Exceeds max number of iterations. Aborted.");
flow_done:
// Calculate three-phase short circuit.
if (!short_circuit)
return;
const auto impedance = calc->node_impedance();
writer->to_csv_file("node-impedance-real.csv", impedance.first);
writer->to_csv_file("node-impedance-imag.csv", impedance.second);
const auto i_f = calc->short_circuit_current();
writer::print_complex("Three-phase short circuit current: ", i_f);
const auto u_f = calc->short_circuit_voltage();
writer->to_csv_file("short-circuit-voltage.csv", u_f, "Ui(real),Ui(imag)");
const auto i = calc->short_circuit_edge_current();
writer->to_csv_file("short-circuit-edge-current.csv", i, "Iij(real),Iij(imag)");
}
}

8
src/executor.hpp
View File

@ -8,14 +8,20 @@
namespace flow
{
/// Forward declaration.
class factory;
/// Controls the execution of this program.
class executor
{
/// The factory instance.
factory* factory_;
public:
/**
* Default constructor.
*/
explicit executor() = default;
explicit executor();
/// Do execute.
void execute(int argc, char** argv) const;

14
src/writer.cpp
View File

@ -25,7 +25,8 @@ namespace flow
const auto width = csbi.dwSize.X;
#else
winsize win;
ioctl(STDOUT_FILENO, TIOCGWINSZ, &win);
if (ioctl(STDOUT_FILENO, TIOCGWINSZ, &win) == -1)
return 10;
const auto width = win.ws_col;
#endif // _WIN32
return std::floor((width - 7) / 11);
@ -58,7 +59,13 @@ namespace flow
});
}
bool writer::to_csv_file(const std::string& path, const arma::mat& mat, const std::string& header)
void writer::print_complex(const std::string& prefix, const std::complex<double>& complex)
{
std::cout << prefix << complex.real() << (complex.imag() < 0 ? '-' : '+') <<
'j' << std::abs(complex.imag()) << std::endl;
}
void writer::to_csv_file(const std::string& path, const arma::mat& mat, const std::string& header) const
{
std::ofstream ofstream;
ofstream.exceptions(std::ifstream::failbit);
@ -83,10 +90,9 @@ namespace flow
}
ofstream << std::endl;
});
return true;
}
catch (const std::exception&) {
return false;
error("Failed to write to file.");
}
}
}

7
src/writer.hpp
View File

@ -83,6 +83,11 @@ namespace flow
*/
static void print_mat(const arma::mat& mat);
/**
* Print a complex number to stdout.
*/
static void print_complex(const std::string& prefix, const std::complex<double>& complex);
/**
* Write a matrix to a file in CSV format.
*
@ -90,6 +95,6 @@ namespace flow
* @param mat Matrix to be printed.
* @param header Header of CSV file
*/
bool to_csv_file(const std::string& path, const arma::mat& mat, const std::string& header = "");
void to_csv_file(const std::string& path, const arma::mat& mat, const std::string& header = "") const;
};
}