Finished most functions.
This commit is contained in:
parent
997daa3426
commit
81ea918884
21
src/args.cpp
21
src/args.cpp
|
@ -1,3 +1,9 @@
|
|||
//
|
||||
// arma-flow/args.cpp
|
||||
//
|
||||
// @author CismonX
|
||||
//
|
||||
|
||||
#include "args.hpp"
|
||||
|
||||
namespace flow
|
||||
|
@ -39,14 +45,21 @@ namespace flow
|
|||
return arg_parser_.getFlag("r");
|
||||
}
|
||||
|
||||
int args::max_iterations()
|
||||
bool args::max_iterations(unsigned& max)
|
||||
{
|
||||
return arg_parser_.getInt("i");
|
||||
const auto arg_i = arg_parser_.getInt("i");
|
||||
if (!arg_parser_.found("i") || arg_i <= 0) {
|
||||
max = 100;
|
||||
return false;
|
||||
}
|
||||
max = arg_i;
|
||||
return true;
|
||||
}
|
||||
|
||||
double args::accuracy()
|
||||
bool args::accuracy(double& epsilon)
|
||||
{
|
||||
return arg_parser_.getDouble("a");
|
||||
epsilon = arg_parser_.getDouble("a");
|
||||
return arg_parser_.found("a");
|
||||
}
|
||||
|
||||
bool args::verbose()
|
||||
|
|
61
src/args.hpp
61
src/args.hpp
|
@ -1,30 +1,85 @@
|
|||
//
|
||||
// arma-flow/args.hpp
|
||||
//
|
||||
// @author CismonX
|
||||
//
|
||||
|
||||
#pragma once
|
||||
|
||||
#include <options.h>
|
||||
|
||||
namespace flow
|
||||
{
|
||||
/// Provides CLI argument parsing utility.
|
||||
class args
|
||||
{
|
||||
/// The argument parser.
|
||||
options::ArgParser arg_parser_;
|
||||
|
||||
public:
|
||||
/**
|
||||
* Constructor.
|
||||
*/
|
||||
explicit args();
|
||||
|
||||
/**
|
||||
* Parse arguments.
|
||||
*
|
||||
* @param argc Number of arguments.
|
||||
* @param argv Arguments.
|
||||
*/
|
||||
void parse(int argc, char** argv);
|
||||
|
||||
/**
|
||||
* Get input file path from arguments.
|
||||
*
|
||||
* @param nodes Path to node data file.
|
||||
* @param edges Path to edge data file.
|
||||
* @return Whether argument is provided.
|
||||
*/
|
||||
bool input_file_path(std::string& nodes, std::string& edges);
|
||||
|
||||
|
||||
/**
|
||||
* Get output file path prefix from arguments.
|
||||
*
|
||||
* @param output Path prefix to output files.
|
||||
* @return Whether argument is provided.
|
||||
*/
|
||||
bool output_file_path(std::string& output);
|
||||
|
||||
/**
|
||||
* Check whether to remove the first line from input files.
|
||||
*
|
||||
* @return Whether argument is provided.
|
||||
*/
|
||||
bool remove_first_line();
|
||||
|
||||
int max_iterations();
|
||||
/**
|
||||
* Check max number of iterations before aborting calculation.
|
||||
*
|
||||
* @return Whether argument is provided.
|
||||
*/
|
||||
bool max_iterations(unsigned& max);
|
||||
|
||||
double accuracy();
|
||||
/**
|
||||
* Check max deviation to be tolerated.
|
||||
*
|
||||
* @return Whether argument is provided.
|
||||
*/
|
||||
bool accuracy(double& epsilon);
|
||||
|
||||
/**
|
||||
* Check whether to enable verbose output.
|
||||
*
|
||||
* @return Whether argument is provided.
|
||||
*/
|
||||
bool verbose();
|
||||
|
||||
/**
|
||||
* Pring help message and exit.
|
||||
*
|
||||
* @return Whether argument is provided.
|
||||
*/
|
||||
void help();
|
||||
};
|
||||
}
|
331
src/calc.cpp
331
src/calc.cpp
|
@ -1,84 +1,323 @@
|
|||
//
|
||||
// arma-flow/calc.cpp
|
||||
//
|
||||
// @author CismonX
|
||||
//
|
||||
|
||||
#include "calc.hpp"
|
||||
#include "executor.hpp"
|
||||
#include "writer.hpp"
|
||||
|
||||
namespace flow
|
||||
{
|
||||
std::complex<double> calc::edge_data::impedance() const
|
||||
unsigned calc::node_offset(unsigned id) const
|
||||
{
|
||||
const auto deno = r * r + x * x;
|
||||
return { r / deno, -x / deno };
|
||||
for (auto offset = 0U; offset < num_nodes_; ++offset)
|
||||
if (nodes_[offset].id == id)
|
||||
return offset;
|
||||
return 0;
|
||||
}
|
||||
|
||||
std::complex<double> calc::edge_data::grounding_admittance() const
|
||||
double calc::j_elem_g_b(unsigned row, unsigned col) const
|
||||
{
|
||||
return { 0, b };
|
||||
return n_adm_g_.at(row, col) * e_[row] + n_adm_b_.at(row, col) * f_[row];
|
||||
}
|
||||
|
||||
double calc::j_elem_b_g(unsigned row, unsigned col) const
|
||||
{
|
||||
return n_adm_b_.at(row, col) * e_[row] - n_adm_g_.at(row, col) * f_[row];
|
||||
}
|
||||
|
||||
double calc::j_elem_a(unsigned row) const
|
||||
{
|
||||
auto sum = 0.0;
|
||||
for (auto col = 0U; col < num_nodes_; ++col)
|
||||
if (row == col || adj_.at(row, col))
|
||||
sum += n_adm_b_.at(row, col) * f_[col] - n_adm_g_.at(row, col) * e_[col];
|
||||
return sum;
|
||||
}
|
||||
|
||||
double calc::j_elem_c(unsigned row) const
|
||||
{
|
||||
auto sum = 0.0;
|
||||
for (auto col = 0U; col < num_nodes_; ++col)
|
||||
if (row == col || adj_.at(row, col))
|
||||
sum += n_adm_g_.at(row, col) * f_[col] + n_adm_b_.at(row, col) * e_[col];
|
||||
return sum;
|
||||
}
|
||||
|
||||
void calc::jacobian()
|
||||
{
|
||||
mat_elem_foreach(j_h_, [this](auto& elem, auto row, auto col)
|
||||
{
|
||||
if (row == col)
|
||||
elem = j_elem_c(row) - j_elem_d(row);
|
||||
else
|
||||
elem = -j_elem_b_g(row, col);
|
||||
});
|
||||
mat_elem_foreach(j_n_, [this](auto& elem, auto row, auto col)
|
||||
{
|
||||
if (row == col)
|
||||
elem = -j_elem_a(row) + j_elem_b(row);
|
||||
else
|
||||
elem = j_elem_g_b(row, col);
|
||||
});
|
||||
mat_elem_foreach(j_m_, [this](auto& elem, auto row, auto col)
|
||||
{
|
||||
if (row == col)
|
||||
elem = -j_elem_a(row) - j_elem_b(row);
|
||||
else
|
||||
elem = -j_elem_g_b(row, col);
|
||||
});
|
||||
mat_elem_foreach(j_l_, [this](auto& elem, auto row, auto col)
|
||||
{
|
||||
if (row == col)
|
||||
elem = -j_elem_c(row) - j_elem_d(row);
|
||||
else
|
||||
elem = -j_elem_b_g(row, col);
|
||||
});
|
||||
// We shouldn't start at 0. We should start at m.
|
||||
mat_elem_foreach(j_r_, [this](auto& elem, auto row, auto col)
|
||||
{
|
||||
if (row + num_pq_ == col)
|
||||
elem = f_[row + num_pq_] * 2;
|
||||
else
|
||||
elem = 0;
|
||||
});
|
||||
mat_elem_foreach(j_s_, [this](auto& elem, auto row, auto col)
|
||||
{
|
||||
if (row + num_pq_ == col)
|
||||
elem = e_[row + num_pq_] * 2;
|
||||
else
|
||||
elem = 0;
|
||||
});
|
||||
}
|
||||
|
||||
void calc::init(const arma::dmat& nodes, const arma::dmat& edges,
|
||||
int max_iterations, double accuracy)
|
||||
bool verbose, double epsilon)
|
||||
{
|
||||
if (max_iterations <= 0)
|
||||
executor::print_and_exit("Bad number of iterations. Positive integer expected.");
|
||||
if (accuracy < 0 || accuracy > 1)
|
||||
executor::print_and_exit("Invalid accuracy.");
|
||||
if (nodes.n_cols != 5 || edges.n_cols != 6)
|
||||
executor::print_and_exit("Bad input matrix format.");
|
||||
writer::error("Bad input matrix format.");
|
||||
nodes.each_row([this](const arma::rowvec& row)
|
||||
{
|
||||
using node_type = node_data::node_type;
|
||||
auto type = node_type::swing;
|
||||
switch (static_cast<unsigned>(row[4]))
|
||||
{
|
||||
case 0:
|
||||
break;
|
||||
case 1:
|
||||
type = node_type::pq;
|
||||
break;
|
||||
case 2:
|
||||
type = node_type::pv;
|
||||
break;
|
||||
default:
|
||||
executor::print_and_exit("Bad node type.");
|
||||
}
|
||||
auto type_val = static_cast<unsigned>(row[4]);
|
||||
auto type = node_data::swing;
|
||||
if (type_val == 1) {
|
||||
type = node_data::pq;
|
||||
++num_pq_;
|
||||
} else if (type_val == 2) {
|
||||
type = node_data::pv;
|
||||
++num_pv_;
|
||||
} else if (type_val != 0)
|
||||
writer::error("Bad node type.");
|
||||
nodes_.push_back({
|
||||
node_num_++, row[0], row[1], row[2], row[3], type
|
||||
num_nodes_++, row[0], row[1], row[2], row[3], type
|
||||
});
|
||||
});
|
||||
// Nodes should be sorted, PQ nodes should be followed by PV nodes,
|
||||
// while swing node be the last.
|
||||
std::sort(nodes_.begin(), nodes_.end(), [](auto&& n1, auto&& n2)
|
||||
{
|
||||
return n1.type < n2.type;
|
||||
});
|
||||
adj_.zeros(num_nodes_, num_nodes_);
|
||||
if (num_nodes_ != num_pq_ + num_pv_ + 1)
|
||||
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;
|
||||
if (n1 >= node_num_ || n2 >= node_num_)
|
||||
executor::print_and_exit("Bad node offset.");
|
||||
if (n1 >= num_nodes_ || n2 >= num_nodes_)
|
||||
writer::error("Bad node offset.");
|
||||
edges_.push_back({
|
||||
&nodes_[n1], &nodes_[n2], row[2], row[3], row[4], row[5]
|
||||
n1, n2, row[2], row[3], row[4], row[5]
|
||||
});
|
||||
});
|
||||
verbose_ = verbose;
|
||||
epsilon_ = epsilon;
|
||||
}
|
||||
|
||||
void calc::node_admittance()
|
||||
{
|
||||
arma::cx_mat node_adm_cplx(node_num_, node_num_, arma::fill::zeros);
|
||||
arma::cx_mat node_adm_cplx(num_nodes_, num_nodes_, arma::fill::zeros);
|
||||
for (auto&& edge : edges_) {
|
||||
const auto n1 = edge.n1->offset;
|
||||
const auto n2 = edge.n2->offset;
|
||||
const auto impedance = edge.impedance();
|
||||
const auto m = node_offset(edge.m);
|
||||
const auto n = node_offset(edge.n);
|
||||
// Whether this edge has transformer.
|
||||
if (edge.k) {
|
||||
node_adm_cplx.at(n1, n1) += impedance;
|
||||
node_adm_cplx.at(n2, n2) += impedance / (edge.k * edge.k);
|
||||
node_adm_cplx.at(n1, n2) -= impedance / edge.k;
|
||||
node_adm_cplx.at(n2, n1) -= impedance / edge.k;
|
||||
}
|
||||
else {
|
||||
node_adm_cplx.at(m, m) += impedance;
|
||||
node_adm_cplx.at(n, n) += impedance / std::pow(edge.k, 2);
|
||||
node_adm_cplx.at(m, n) -= impedance / edge.k;
|
||||
node_adm_cplx.at(n, m) -= impedance / edge.k;
|
||||
} else {
|
||||
const auto delta_diag = impedance + edge.grounding_admittance();
|
||||
node_adm_cplx.at(n1, n1) += delta_diag;
|
||||
node_adm_cplx.at(n2, n2) += delta_diag;
|
||||
node_adm_cplx.at(n1, n2) -= impedance;
|
||||
node_adm_cplx.at(n2, n1) -= impedance;
|
||||
node_adm_cplx.at(m, m) += delta_diag;
|
||||
node_adm_cplx.at(n, n) += delta_diag;
|
||||
node_adm_cplx.at(m, n) -= impedance;
|
||||
node_adm_cplx.at(n, m) -= impedance;
|
||||
}
|
||||
adj_.at(m, n) = 1;
|
||||
adj_.at(n, m) = 1;
|
||||
}
|
||||
node_adm_real_ = arma::real(node_adm_cplx);
|
||||
node_adm_imag_ = arma::imag(node_adm_cplx);
|
||||
n_adm_g_ = arma::real(node_adm_cplx);
|
||||
n_adm_b_ = arma::imag(node_adm_cplx);
|
||||
if (verbose_) {
|
||||
writer::println("Real part of node admittance matrix:");
|
||||
writer::print_mat(n_adm_g_);
|
||||
writer::println("Imaginary part of node admittance matrix:");
|
||||
writer::print_mat(n_adm_b_);
|
||||
}
|
||||
}
|
||||
|
||||
void calc::iterate_init()
|
||||
{
|
||||
init_p_.zeros(num_nodes_ - 1);
|
||||
init_q_.zeros(num_pq_);
|
||||
init_v_.zeros(num_pv_);
|
||||
auto i_p = 0;
|
||||
auto i_q = 0;
|
||||
auto i_v = 0;
|
||||
e_.resize(num_nodes_);
|
||||
for (auto&& node : nodes_) {
|
||||
if (node.type == node_data::pq) {
|
||||
init_p_[i_p] = -node.p;
|
||||
init_q_[i_q++] = -node.q;
|
||||
} else if (node.type == node_data::pv) {
|
||||
init_p_[i_p] = node.g - node.p;
|
||||
init_v_[i_v++] = node.v;
|
||||
}
|
||||
e_[i_p++] = node.v;
|
||||
}
|
||||
f_.zeros(num_nodes_);
|
||||
j_h_.zeros(num_nodes_ - 1, num_nodes_ - 1);
|
||||
j_n_.zeros(num_nodes_ - 1, num_nodes_ - 1);
|
||||
j_m_.zeros(num_pq_, num_nodes_ - 1);
|
||||
j_l_.zeros(num_pq_, num_nodes_ - 1);
|
||||
j_r_.zeros(num_pv_, num_nodes_ - 1);
|
||||
j_s_.zeros(num_pv_, num_nodes_ - 1);
|
||||
delta_p_.zeros(num_nodes_ - 1);
|
||||
delta_q_.zeros(num_pq_);
|
||||
delta_v_.zeros(num_pv_);
|
||||
p_.zeros(num_nodes_);
|
||||
q_.zeros(num_nodes_);
|
||||
update_f_x();
|
||||
}
|
||||
|
||||
void calc::prepare_solve()
|
||||
{
|
||||
// Cross-construct F(x) vector.
|
||||
f_x_.zeros(2 * num_nodes_ - 2);
|
||||
for (auto row = 0U; row < num_nodes_ - 1; ++row) {
|
||||
f_x_[2 * row] = delta_p_[row];
|
||||
f_x_[2 * row + 1] = row < num_pq_ ? delta_q_[row] : delta_v_[row - num_pq_];
|
||||
}
|
||||
// Cross-construct jacobian matrix.
|
||||
arma::dmat jacobian(2 * num_nodes_ - 2, 2 * num_nodes_ - 2);
|
||||
for (auto row = 0U; row < num_nodes_ - 1; ++row)
|
||||
for (auto col = 0U; col < num_nodes_ - 1; ++col) {
|
||||
jacobian.at(2 * row, 2 * col) = j_h_.at(row, col);
|
||||
jacobian.at(2 * row, 2 * col + 1) = j_n_.at(row, col);
|
||||
if (row < num_pq_) {
|
||||
jacobian.at(2 * row + 1, 2 * col) = j_m_.at(row, col);
|
||||
jacobian.at(2 * row + 1, 2 * col + 1) = j_l_.at(row, col);
|
||||
} else {
|
||||
jacobian.at(2 * row + 1, 2 * col) = j_r_.at(row - num_pq_, col);
|
||||
jacobian.at(2 * row + 1, 2 * col + 1) = j_s_.at(row - num_pq_, col);
|
||||
}
|
||||
}
|
||||
if (verbose_) {
|
||||
writer::println("Jacobian matrix");
|
||||
writer::print_mat(jacobian);
|
||||
}
|
||||
// Transform into a sparse matrix for spsolve().
|
||||
j_ = jacobian;
|
||||
}
|
||||
|
||||
void calc::update_f_x()
|
||||
{
|
||||
vec_elem_foreach(delta_p_, [this](auto& elem, auto row)
|
||||
{
|
||||
p_[row] = calc_p(row);
|
||||
elem = init_p_[row] - p_[row];
|
||||
});
|
||||
vec_elem_foreach(delta_q_, [this](auto& elem, auto row)
|
||||
{
|
||||
q_[row] = calc_q(row);
|
||||
elem = init_q_[row] - q_[row];
|
||||
});
|
||||
vec_elem_foreach(delta_v_, [this](auto& elem, auto row)
|
||||
{
|
||||
auto i = row + num_pq_;
|
||||
elem = std::pow(init_v_[row], 2) - std::pow(e_[i], 2) - std::pow(f_[i], 2);
|
||||
});
|
||||
if (verbose_) {
|
||||
writer::println("Delta P:");
|
||||
writer::print_mat(delta_p_.t());
|
||||
writer::println("Delta Q:");
|
||||
writer::print_mat(delta_q_.t());
|
||||
writer::println("Delta U^2:");
|
||||
writer::print_mat(delta_v_.t());
|
||||
}
|
||||
}
|
||||
|
||||
unsigned calc::solve()
|
||||
{
|
||||
jacobian();
|
||||
prepare_solve();
|
||||
const auto x_vec = spsolve(j_, f_x_, "superlu").col(0);
|
||||
vec_elem_foreach(f_, [&x_vec](auto& elem, auto row)
|
||||
{
|
||||
elem += x_vec[2 * row];
|
||||
});
|
||||
vec_elem_foreach(e_, [&x_vec](auto& elem, auto row)
|
||||
{
|
||||
elem += x_vec[2 * row + 1];
|
||||
});
|
||||
if (verbose_) {
|
||||
writer::println("Correction vector of voltage (real):");
|
||||
writer::print_mat(e_.t());
|
||||
writer::println("Correction vector of voltage (imaginary):");
|
||||
writer::print_mat(f_.t());
|
||||
}
|
||||
update_f_x();
|
||||
writer::println("Number of iterations: ", ++n_iter_);
|
||||
return n_iter_;
|
||||
}
|
||||
|
||||
double calc::get_max() const
|
||||
{
|
||||
arma::dmat mat = join_cols(join_cols(delta_p_, delta_q_), delta_v_);
|
||||
auto max = 0.0;
|
||||
for (auto&& elem : mat)
|
||||
if (std::abs(elem) > max)
|
||||
max = std::abs(elem);
|
||||
return max;
|
||||
}
|
||||
|
||||
arma::dmat calc::result()
|
||||
{
|
||||
p_[num_nodes_ - 1] = calc_p(num_nodes_ - 1);
|
||||
vec_elem_foreach(delta_v_, [this](auto& elem, auto row)
|
||||
{
|
||||
auto i = row + num_pq_;
|
||||
q_[i] = calc_q(i);
|
||||
});
|
||||
q_[num_nodes_ - 1] = calc_q(num_nodes_ - 1);
|
||||
arma::colvec v(num_nodes_);
|
||||
vec_elem_foreach(v, [this](auto& elem, auto col)
|
||||
{
|
||||
elem = std::sqrt(std::pow(e_[col], 2) + std::pow(f_[col], 2));
|
||||
});
|
||||
arma::colvec theta(num_nodes_);
|
||||
vec_elem_foreach(theta, [this](auto& elem, auto col)
|
||||
{
|
||||
elem = std::atan(f_[col] / e_[col]);
|
||||
});
|
||||
arma::colvec node_id(num_nodes_);
|
||||
vec_elem_foreach(node_id, [this](auto& elem, auto col)
|
||||
{
|
||||
elem = nodes_[col].id;
|
||||
});
|
||||
return join_rows(node_id, join_rows(join_rows(p_, q_), join_rows(v, theta)));
|
||||
}
|
||||
}
|
||||
|
|
238
src/calc.hpp
238
src/calc.hpp
|
@ -1,3 +1,9 @@
|
|||
//
|
||||
// arma-flow/calc.hpp
|
||||
//
|
||||
// @author CismonX
|
||||
//
|
||||
|
||||
#pragma once
|
||||
|
||||
#include <armadillo>
|
||||
|
@ -5,51 +11,245 @@
|
|||
|
||||
namespace flow
|
||||
{
|
||||
/// Power flow calculation.
|
||||
class calc
|
||||
{
|
||||
/// Structure of node data.
|
||||
struct node_data
|
||||
{
|
||||
unsigned long offset;
|
||||
double v, g, p, q;
|
||||
/// Node ID.
|
||||
unsigned id;
|
||||
|
||||
/// Voltage (real);
|
||||
double v;
|
||||
|
||||
/// Generator (active power)
|
||||
double g;
|
||||
|
||||
/// Load (active power)
|
||||
double p;
|
||||
|
||||
/// Load (reactive power)
|
||||
double q;
|
||||
|
||||
/// Node type.
|
||||
enum node_type {
|
||||
pq, pv, swing
|
||||
} type;
|
||||
};
|
||||
|
||||
/// Structure of edge data.
|
||||
struct edge_data
|
||||
{
|
||||
node_data *n1, *n2;
|
||||
double r, x, b, k;
|
||||
std::complex<double> impedance() const;
|
||||
std::complex<double> grounding_admittance() const;
|
||||
/// First and second node ID.
|
||||
unsigned m, n;
|
||||
|
||||
/// Resistance (real).
|
||||
double r;
|
||||
|
||||
/// Resistance (imaginary).
|
||||
double x;
|
||||
|
||||
/// Gounding admittance (imaginary).
|
||||
double b;
|
||||
|
||||
/// Transformer ratio.
|
||||
double k;
|
||||
|
||||
/**
|
||||
* Get impedance of edge.
|
||||
*
|
||||
* @return Impedance (complex).
|
||||
*/
|
||||
std::complex<double> impedance() const
|
||||
{
|
||||
const auto deno = r * r + x * x;
|
||||
return { r / deno, -x / deno };
|
||||
}
|
||||
|
||||
/**
|
||||
* Get gounding admittance.
|
||||
*
|
||||
* @return Grounding admittance (complex).
|
||||
*/
|
||||
std::complex<double> grounding_admittance() const
|
||||
{
|
||||
return { 0, b };
|
||||
};
|
||||
};
|
||||
|
||||
/// Vector of nodes.
|
||||
std::vector<node_data> nodes_;
|
||||
|
||||
unsigned long node_num_ = 0;
|
||||
/// Number of nodes.
|
||||
unsigned num_nodes_ = 0;
|
||||
|
||||
/// Number of PQ nodes and PV nodes.
|
||||
unsigned num_pq_ = 0, num_pv_ = 0;
|
||||
|
||||
/// Vector of edges.
|
||||
std::vector<edge_data> edges_;
|
||||
|
||||
arma::dmat node_adm_real_;
|
||||
/// Adjacency matrix of nodes.
|
||||
arma::uchar_mat adj_;
|
||||
|
||||
arma::dmat node_adm_imag_;
|
||||
/// Node admittance matrix.
|
||||
arma::dmat n_adm_g_, n_adm_b_;
|
||||
|
||||
/// Given values of power and voltage.
|
||||
arma::colvec init_p_, init_q_, init_v_;
|
||||
|
||||
/// Correction vector of voltage.
|
||||
arma::colvec e_, f_;
|
||||
|
||||
/// Imbalace of active/reactive power and voltage.
|
||||
arma::colvec delta_p_, delta_q_, delta_v_;
|
||||
|
||||
/// Submatrix of Jacobian matrix.
|
||||
arma::dmat j_h_, j_n_, j_m_, j_l_, j_r_, j_s_;
|
||||
|
||||
/// F(x) of jacobian matrix.
|
||||
arma::colvec f_x_;
|
||||
|
||||
/// Jacobian matrix (sparse).
|
||||
arma::sp_dmat j_;
|
||||
|
||||
/// Power vector of nodes.
|
||||
arma::colvec p_, q_;
|
||||
|
||||
/// Whether verbose output is enabled.
|
||||
bool verbose_ = false;
|
||||
|
||||
/// Max deviation to be tolerated.
|
||||
double epsilon_;
|
||||
|
||||
/// Number of iterations.
|
||||
unsigned n_iter_ = 0;
|
||||
|
||||
/**
|
||||
* Get offset of sorted node by ID.
|
||||
*
|
||||
* @param id Node ID.
|
||||
* @return Node offset.
|
||||
*/
|
||||
unsigned node_offset(unsigned id) const;
|
||||
|
||||
/// G(i, j) * e(i) + B(i, j) * f(i)
|
||||
double j_elem_g_b(unsigned row, unsigned col) const;
|
||||
|
||||
/// B(i, j) * e(i) - G(i, j) * f(i)
|
||||
double j_elem_b_g(unsigned row, unsigned col) const;
|
||||
|
||||
/// G(i, i) * e(i) + B(i, i) * f(i)
|
||||
double j_elem_b(unsigned row) const
|
||||
{
|
||||
return j_elem_g_b(row, row);
|
||||
}
|
||||
|
||||
/// B(i, i) * e(i) - G(i, i) * f(i)
|
||||
double j_elem_d(unsigned row) const
|
||||
{
|
||||
return j_elem_b_g(row, row);
|
||||
}
|
||||
|
||||
/// sum(i == j || adj(i, j), B(i, j) * f(j) - G(i, j) * e(j))
|
||||
double j_elem_a(unsigned row) const;
|
||||
|
||||
/// sum(i == j || adj(i, j), G(i, j) * f(j) + B(i, j) * e(j))
|
||||
double j_elem_c(unsigned row) const;
|
||||
|
||||
/**
|
||||
* Traverse a matrix.
|
||||
*
|
||||
* @param mat Matrix to be traversed.
|
||||
* @param func Callback for each element.
|
||||
*/
|
||||
template <typename M, typename F>
|
||||
static void mat_elem_foreach(M& mat, F func)
|
||||
{
|
||||
for (auto row = 0U; row < mat.n_rows; ++row)
|
||||
for (auto col = 0U; col < mat.n_cols; ++col)
|
||||
func(mat.at(row, col), row, col);
|
||||
}
|
||||
|
||||
/**
|
||||
* Traverse a vector.
|
||||
*
|
||||
* @param vec Vector to be traversed.
|
||||
* @param func Callback for each element
|
||||
*/
|
||||
template <typename V, typename F>
|
||||
static void vec_elem_foreach(V& vec, F func)
|
||||
{
|
||||
for (auto i = 0U; i < vec.n_elem; ++i)
|
||||
func(vec[i], i);
|
||||
}
|
||||
|
||||
/**
|
||||
* Update F(x) of jacobian matrix.
|
||||
*/
|
||||
void update_f_x();
|
||||
|
||||
/**
|
||||
* Calculate jacobian matrix.
|
||||
*/
|
||||
void jacobian();
|
||||
|
||||
/**
|
||||
* Prepare to solve the formula.
|
||||
*/
|
||||
void prepare_solve();
|
||||
|
||||
/// Calculate active power of a node.
|
||||
double calc_p(unsigned row) const
|
||||
{
|
||||
return f_[row] * j_elem_c(row) - e_[row] * j_elem_a(row);
|
||||
}
|
||||
|
||||
/// Calculate reactive power of a node.
|
||||
double calc_q(unsigned row) const
|
||||
{
|
||||
return -f_[row] * j_elem_a(row) - e_[row] * j_elem_c(row);
|
||||
}
|
||||
|
||||
public:
|
||||
/**
|
||||
* Default constructor.
|
||||
*/
|
||||
explicit calc() = default;
|
||||
|
||||
/**
|
||||
* Initialize.
|
||||
*/
|
||||
void init(const arma::dmat& nodes, const arma::dmat& edges,
|
||||
int max_iterations, double accuracy);
|
||||
bool verbose, double epsilon);
|
||||
|
||||
/**
|
||||
* Calculate node admittance.
|
||||
*/
|
||||
void node_admittance();
|
||||
|
||||
arma::dmat* get_node_adm_real()
|
||||
{
|
||||
return &node_adm_real_;
|
||||
}
|
||||
/**
|
||||
* Initialize iteratiob.
|
||||
*/
|
||||
void iterate_init();
|
||||
|
||||
arma::dmat* get_node_adm_imag()
|
||||
{
|
||||
return &node_adm_imag_;
|
||||
}
|
||||
/**
|
||||
* Solve formula.
|
||||
*
|
||||
* @return number of iterations.
|
||||
*/
|
||||
unsigned solve();
|
||||
|
||||
/**
|
||||
* Get max value of correction vector.
|
||||
*/
|
||||
double get_max() const;
|
||||
|
||||
/**
|
||||
* Get result of calculation
|
||||
*
|
||||
*/
|
||||
arma::dmat result();
|
||||
};
|
||||
}
|
||||
}
|
||||
|
|
|
@ -1,49 +1,73 @@
|
|||
//
|
||||
// arma-flow/executor.cpp
|
||||
//
|
||||
// @author CismonX
|
||||
//
|
||||
|
||||
#include <iostream>
|
||||
|
||||
#include "executor.hpp"
|
||||
#include "factory.hpp"
|
||||
#include "output.hpp"
|
||||
#include "writer.hpp"
|
||||
|
||||
namespace flow
|
||||
{
|
||||
void executor::print_and_exit(const std::string& message)
|
||||
{
|
||||
std::cout << message << std::endl;
|
||||
exit(1);
|
||||
}
|
||||
|
||||
void executor::execute(int argc, char** argv) const
|
||||
{
|
||||
// Get componenets.
|
||||
auto factory = factory::get();
|
||||
auto args = factory->get_args();
|
||||
auto input = factory->get_reader();
|
||||
auto calc = factory->get_calc();
|
||||
|
||||
// Parse options.
|
||||
args->parse(argc, argv);
|
||||
|
||||
// Read data from file.
|
||||
std::string path_to_nodes, path_to_edges;
|
||||
if (!args->input_file_path(path_to_nodes, path_to_edges))
|
||||
args->help();
|
||||
auto input = factory->get_input();
|
||||
const auto remove = args->remove_first_line();
|
||||
if (!input->from_csv_file(path_to_nodes, remove))
|
||||
print_and_exit("Failed to read edge data from file.");
|
||||
writer::error("Failed to read node data from file.");
|
||||
const auto nodes = input->get_mat();
|
||||
if (!input->from_csv_file(path_to_edges, remove))
|
||||
print_and_exit("Failed to read edge data from file.");
|
||||
writer::error("Failed to read edge data from file.");
|
||||
const auto edges = input->get_mat();
|
||||
auto calc = factory->get_calc();
|
||||
calc->init(nodes, edges, args->max_iterations(), args->accuracy());
|
||||
calc->node_admittance();
|
||||
const auto node_adm_real = calc->get_node_adm_real();
|
||||
const auto node_adm_imag = calc->get_node_adm_imag();
|
||||
const auto verbose = args->verbose();
|
||||
if (verbose) {
|
||||
std::cout << "Real part of node admittance matrix:" << std::endl;
|
||||
output::print_dmat(*node_adm_real);
|
||||
std::cout << "Imaginary part of node admittance matrix:" << std::endl;
|
||||
output::print_dmat(*node_adm_imag);
|
||||
}
|
||||
std::string output_path;
|
||||
if (!args->output_file_path(output_path))
|
||||
std::cout << "Output file path not specified. Using result-*.csv by default." << std::endl;
|
||||
|
||||
// Get options.
|
||||
const auto verbose = args->verbose();
|
||||
auto max = 0U;
|
||||
if (!args->max_iterations(max) && verbose)
|
||||
writer::notice("Max number of iterations not specified. Defaulted to 100.");
|
||||
auto epsilon = 0.0;
|
||||
if (!args->accuracy(epsilon) && verbose)
|
||||
writer::notice("Accuracy not specified. Defaulted to 0.00001.");
|
||||
if (epsilon < 0 || epsilon > 1)
|
||||
writer::error("Invalid accuracy.");
|
||||
std::string output_path;
|
||||
if (!args->output_file_path(output_path) && verbose)
|
||||
writer::notice("Output file path not specified. Defaulted to result-*.csv.");
|
||||
|
||||
// Initialize calculation.
|
||||
calc->init(nodes, edges, verbose, epsilon);
|
||||
calc->node_admittance();
|
||||
calc->iterate_init();
|
||||
|
||||
// Do iteration.
|
||||
unsigned num_iterations;
|
||||
do {
|
||||
num_iterations = calc->solve();
|
||||
if (calc->get_max() <= epsilon) {
|
||||
writer::println("Finished. Total number of iterations: ", num_iterations);
|
||||
const auto result = calc->result();
|
||||
if (verbose) {
|
||||
writer::println("Result(n, P, Q, V, theta):");
|
||||
writer::print_mat(result);
|
||||
}
|
||||
break;
|
||||
}
|
||||
} while (num_iterations < max);
|
||||
writer::println("Exceeds max number of iterations. Aborted.");
|
||||
}
|
||||
}
|
||||
|
|
|
@ -1,15 +1,22 @@
|
|||
#pragma once
|
||||
//
|
||||
// arma-flow/executor.hpp
|
||||
//
|
||||
// @author CismonX
|
||||
//
|
||||
|
||||
#pragma once
|
||||
|
||||
namespace flow
|
||||
{
|
||||
/// Controls the excution of this program.
|
||||
class executor
|
||||
{
|
||||
public:
|
||||
/**
|
||||
* Default constructor.
|
||||
*/
|
||||
explicit executor() = default;
|
||||
|
||||
static void print_and_exit(const std::string& message);
|
||||
|
||||
/// Do execute.
|
||||
void execute(int argc, char** argv) const;
|
||||
};
|
||||
}
|
|
@ -1,61 +1,82 @@
|
|||
//
|
||||
// arma-flow/factory.hpp
|
||||
//
|
||||
// @author CismonX
|
||||
//
|
||||
|
||||
#pragma once
|
||||
|
||||
#include "input.hpp"
|
||||
#include "reader.hpp"
|
||||
#include "args.hpp"
|
||||
#include "executor.hpp"
|
||||
#include "calc.hpp"
|
||||
#include "output.hpp"
|
||||
#include "writer.hpp"
|
||||
|
||||
namespace flow
|
||||
{
|
||||
/// Contains factory methods for this program.
|
||||
class factory
|
||||
{
|
||||
input input_;
|
||||
/// The reader instance.
|
||||
reader reader_;
|
||||
|
||||
/// The arg parser instance.
|
||||
args args_;
|
||||
|
||||
/// The executor instance.
|
||||
executor executor_;
|
||||
|
||||
/// The power flow calculator instance.
|
||||
calc calc_;
|
||||
|
||||
output output_;
|
||||
/// The writer instance.
|
||||
writer writer_;
|
||||
|
||||
/// Singleton for this factory.
|
||||
static factory singleton_;
|
||||
|
||||
/**
|
||||
* Private default constructor.
|
||||
*/
|
||||
explicit factory() = default;
|
||||
|
||||
public:
|
||||
|
||||
public:
|
||||
/// Get factory.
|
||||
static factory* get()
|
||||
{
|
||||
return &singleton_;
|
||||
}
|
||||
|
||||
input* get_input()
|
||||
/// Get reader.
|
||||
reader* get_reader()
|
||||
{
|
||||
return &input_;
|
||||
return &reader_;
|
||||
}
|
||||
|
||||
/// Get arg parser.
|
||||
args* get_args()
|
||||
{
|
||||
return &args_;
|
||||
}
|
||||
|
||||
/// Get executor.
|
||||
executor* get_executor()
|
||||
{
|
||||
return &executor_;
|
||||
}
|
||||
|
||||
/// Get power flow calculator.
|
||||
calc* get_calc()
|
||||
{
|
||||
return &calc_;
|
||||
}
|
||||
|
||||
output* get_output()
|
||||
/// Get writer.
|
||||
writer* get_writer()
|
||||
{
|
||||
return &output_;
|
||||
return &writer_;
|
||||
}
|
||||
};
|
||||
|
||||
inline factory factory::singleton_;
|
||||
}
|
||||
}
|
||||
|
|
|
@ -1,6 +1,13 @@
|
|||
//
|
||||
// arma-flow/main.cpp
|
||||
//
|
||||
// @author CismonX
|
||||
//
|
||||
|
||||
#include "factory.hpp"
|
||||
|
||||
/// Bootstrap the program.
|
||||
int main(int argc, char** argv)
|
||||
{
|
||||
flow::factory::get()->get_executor()->execute(argc, argv);
|
||||
}
|
||||
}
|
||||
|
|
|
@ -1,16 +0,0 @@
|
|||
#pragma once
|
||||
|
||||
#include <armadillo>
|
||||
|
||||
|
||||
namespace flow
|
||||
{
|
||||
class output
|
||||
{
|
||||
static int max_elems_per_line();
|
||||
public:
|
||||
|
||||
static void print_dmat(const arma::dmat& mat);
|
||||
|
||||
};
|
||||
}
|
|
@ -1,16 +1,28 @@
|
|||
//
|
||||
// arma-flow/reader.hpp
|
||||
//
|
||||
// @author CismonX
|
||||
//
|
||||
|
||||
#include <experimental/filesystem>
|
||||
|
||||
#include "input.hpp"
|
||||
#include "reader.hpp"
|
||||
|
||||
namespace flow
|
||||
{
|
||||
bool input::from_csv_file(const std::string& path, bool remove_first_line)
|
||||
bool reader::from_csv_file(const std::string& path, bool remove_first_line)
|
||||
{
|
||||
std::ifstream ifstream;
|
||||
ifstream.exceptions(std::ifstream::failbit);
|
||||
namespace fs = std::experimental::filesystem;
|
||||
try {
|
||||
ifstream.open(path[0] == '/' ? path : fs::current_path().string() + '/' + path);
|
||||
ifstream.open(
|
||||
#ifdef _WIN32
|
||||
path[1] == ':'
|
||||
#else
|
||||
path[0] == '/'
|
||||
#endif // _WIN32
|
||||
? path : fs::current_path().string() + '/' + path);
|
||||
if (remove_first_line)
|
||||
ifstream.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
|
||||
}
|
||||
|
@ -20,7 +32,7 @@ namespace flow
|
|||
return do_read(ifstream);
|
||||
}
|
||||
|
||||
const arma::dmat& input::get_mat() const
|
||||
const arma::dmat& reader::get_mat() const
|
||||
{
|
||||
return mat_;
|
||||
}
|
|
@ -1,10 +1,17 @@
|
|||
//
|
||||
// arma-flow/reader.hpp
|
||||
//
|
||||
// @author CismonX
|
||||
//
|
||||
|
||||
#pragma once
|
||||
|
||||
#include <armadillo>
|
||||
|
||||
namespace flow
|
||||
{
|
||||
class input
|
||||
/// Provides read utilities.
|
||||
class reader
|
||||
{
|
||||
/// The loaded matrix.
|
||||
arma::dmat mat_;
|
|
@ -1,7 +1,8 @@
|
|||
#include <iomanip>
|
||||
|
||||
#include "output.hpp"
|
||||
#include "executor.hpp"
|
||||
//
|
||||
// arma-flow/writer.cpp
|
||||
//
|
||||
// @author CismonX
|
||||
//
|
||||
|
||||
#ifdef _WIN32
|
||||
#include <Windows.h>
|
||||
|
@ -9,11 +10,14 @@
|
|||
#include <sys/ioctl.h>
|
||||
#include <unistd.h>
|
||||
#endif // _WIN32
|
||||
#include <iomanip>
|
||||
|
||||
#include "writer.hpp"
|
||||
#include "executor.hpp"
|
||||
|
||||
namespace flow
|
||||
{
|
||||
int output::max_elems_per_line()
|
||||
int writer::max_elems_per_line()
|
||||
{
|
||||
#ifdef _WIN32
|
||||
CONSOLE_SCREEN_BUFFER_INFO csbi;
|
||||
|
@ -28,11 +32,11 @@ namespace flow
|
|||
return std::floor((width - 7) / 11);
|
||||
}
|
||||
|
||||
void output::print_dmat(const arma::dmat& mat)
|
||||
void writer::print_mat(const arma::dmat& mat)
|
||||
{
|
||||
mat.each_row([](const arma::rowvec& row)
|
||||
{
|
||||
std::cout << std::setprecision(7) << std::left;
|
||||
std::cout << std::setprecision(6) << std::left;
|
||||
auto counter = 0;
|
||||
auto elems = max_elems_per_line();
|
||||
for (auto&& elem : row)
|
||||
|
@ -41,7 +45,7 @@ namespace flow
|
|||
std::cout << "...(" << row.n_elem - elems << ')';
|
||||
break;
|
||||
}
|
||||
std::cout << std::setw(10) << elem;
|
||||
std::cout << std::setw(10) << elem << ' ';
|
||||
}
|
||||
std::cout << std::endl;
|
||||
});
|
|
@ -0,0 +1,65 @@
|
|||
//
|
||||
// arma-flow/writer.hpp
|
||||
//
|
||||
// @author CismonX
|
||||
//
|
||||
|
||||
#pragma once
|
||||
|
||||
#include <armadillo>
|
||||
|
||||
namespace flow
|
||||
{
|
||||
/// Provides write utilities.
|
||||
class writer
|
||||
{
|
||||
/**
|
||||
* Determines width of stdout.
|
||||
*
|
||||
* @return Character width.
|
||||
*/
|
||||
static int max_elems_per_line();
|
||||
|
||||
public:
|
||||
/**
|
||||
* Print a line to stdout.
|
||||
*
|
||||
* @param message Messages to be printed.
|
||||
*/
|
||||
template <typename ...T>
|
||||
static void println(T&&... message)
|
||||
{
|
||||
(std::cout << ... << message) << std::endl;
|
||||
}
|
||||
|
||||
/**
|
||||
* Print a notice message to stdout.
|
||||
*
|
||||
* @param message Messages to be printed.
|
||||
*/
|
||||
template <typename ...T>
|
||||
static void notice(T&&... message)
|
||||
{
|
||||
(std::cout << "Notice: " << ... << message) << std::endl;
|
||||
}
|
||||
|
||||
/**
|
||||
* Print a error message to stdout and terminate proogram.
|
||||
*
|
||||
* @param message Messages to be printed.
|
||||
*/
|
||||
template <typename ...T>
|
||||
static void error(T&&... message)
|
||||
{
|
||||
(std::cout << "Error: " << ... << message) << std::endl;
|
||||
exit(1);
|
||||
}
|
||||
|
||||
/**
|
||||
* Print a matrix to stdout.
|
||||
*
|
||||
* @param mat Matrix to be printed.
|
||||
*/
|
||||
static void print_mat(const arma::dmat& mat);
|
||||
};
|
||||
}
|
Reference in New Issue