Skip to content

Commit

Permalink
compiles and basic checks work
Browse files Browse the repository at this point in the history
  • Loading branch information
@BDonnot
BDonnot committed Dec 13, 2023
1 parent 8334ef7 commit cc25e53
Show file tree
Hide file tree
Showing 18 changed files with 354 additions and 132 deletions.
5 changes: 2 additions & 3 deletions lightsim2grid/lightSimBackend.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -894,9 +894,9 @@ def runpf(self, is_dc=False):
self._grid.deactivate_result_computation()
# if I init with dc values, it should depends on previous state
self.V[:] = self._grid.get_init_vm_pu() # see issue 30
print("before dc pf")
# print(f"{self.V[:14] = }")
Vdc = self._grid.dc_pf(copy.deepcopy(self.V), self.max_it, self.tol)
print("after dc pf")
# print(f"{Vdc[:14] = }")
self._grid.reactivate_result_computation()
if Vdc.shape[0] == 0:
raise BackendError(f"Divergence of DC powerflow (non connected grid) at the initialization of AC powerflow. Detailed error: {self._grid.get_dc_solver().get_error()}")
Expand All @@ -905,7 +905,6 @@ def runpf(self, is_dc=False):
V_init = copy.deepcopy(self.V)
tick = time.perf_counter()
self._timer_preproc += tick - beg_preproc
print("before ac pf")
V = self._grid.ac_pf(V_init, self.max_it, self.tol)
self._timer_solver += time.perf_counter() - tick
if V.shape[0] == 0:
Expand Down
3 changes: 2 additions & 1 deletion setup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -157,7 +157,8 @@
"src/BaseMultiplePowerflow.cpp",
"src/Computers.cpp",
"src/SecurityAnalysis.cpp",
"src/Solvers.cpp"]
"src/Solvers.cpp",
"src/Utils.cpp"]

if KLU_SOLVER_AVAILABLE:
src_files.append("src/KLUSolver.cpp")
Expand Down
48 changes: 32 additions & 16 deletions src/BaseNRSolver.tpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,15 +13,15 @@

template<class LinearSolver>
bool BaseNRSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type> & Ybus,
CplxVect & V,
const CplxVect & Sbus,
const Eigen::VectorXi & slack_ids,
const RealVect & slack_weights,
const Eigen::VectorXi & pv,
const Eigen::VectorXi & pq,
int max_iter,
real_type tol
)
CplxVect & V,
const CplxVect & Sbus,
const Eigen::VectorXi & slack_ids,
const RealVect & slack_weights,
const Eigen::VectorXi & pv,
const Eigen::VectorXi & pq,
int max_iter,
real_type tol
)
{
/**
This method uses the newton raphson algorithm to compute voltage angles and magnitudes at each bus
Expand All @@ -36,19 +36,28 @@ bool BaseNRSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
// TODO DEBUG MODE
std::ostringstream exc_;
exc_ << "BaseNRSolver::compute_pf: Size of the Sbus should be the same as the size of Ybus. Currently: ";
exc_ << "Sbus (" << Sbus.size() << ") and Ybus (" << Ybus.rows() << ", " << Ybus.rows() << ").";
exc_ << "Sbus (" << Sbus.size() << ") and Ybus (" << Ybus.rows() << ", " << Ybus.cols() << ").";
throw std::runtime_error(exc_.str());
}
if(V.size() != Ybus.rows() || V.size() != Ybus.cols() ){
// TODO DEBUG MODE
std::ostringstream exc_;
exc_ << "BaseNRSolver::compute_pf: Size of V (init voltages) should be the same as the size of Ybus. Currently: ";
exc_ << "V (" << V.size() << ") and Ybus (" << Ybus.rows()<< ", " << Ybus.rows() << ").";
exc_ << "V (" << V.size() << ") and Ybus (" << Ybus.rows()<< ", " << Ybus.cols() << ").";
throw std::runtime_error(exc_.str());
}
reset_timer();
// std::cout << "dist slack" << std::endl;

if(_solver_control.need_reset_solver() ||
_solver_control.has_dimension_changed()){
reset();
}
auto timer = CustTimer();
if(!is_linear_solver_valid()) return false;
if(!is_linear_solver_valid()) {
// err_ = ErrorType::NotInitError;
return false;
}

err_ = ErrorType::NoError; // reset the error if previous error happened

Expand Down Expand Up @@ -85,6 +94,11 @@ bool BaseNRSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
bool has_just_been_initialized = false; // to avoid a call to klu_refactor follow a call to klu_factor in the same loop
// std::cout << "iter " << nr_iter_ << " dx(0): " << -F(0) << " dx(1): " << -F(1) << std::endl;
// std::cout << "slack_absorbed " << slack_absorbed << std::endl;
BaseNRSolver<LinearSolver>::value_map_.clear(); // TODO smarter solver: only needed if ybus has changed
BaseNRSolver<LinearSolver>::dS_dVm_.resize(0,0); // TODO smarter solver: only needed if ybus has changed
BaseNRSolver<LinearSolver>::dS_dVa_.resize(0,0); // TODO smarter solver: only needed if ybus has changed
// BaseNRSolver<LinearSolver>::dS_dVm_.setZero(); // TODO smarter solver: only needed if ybus has changed
// BaseNRSolver<LinearSolver>::dS_dVa_.setZero(); // TODO smarter solver: only needed if ybus has changed
while ((!converged) & (nr_iter_ < max_iter)){
nr_iter_++;
fill_jacobian_matrix(Ybus, V_, slack_bus_id, slack_weights, pq, pvpq, pq_inv, pvpq_inv);
Expand Down Expand Up @@ -146,6 +160,7 @@ bool BaseNRSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
<< "\n\t timer_total_nr_: " << timer_total_nr_
<< "\n\n";
#endif // __COUT_TIMES
_solver_control.tell_none_changed();
return res;
}

Expand All @@ -165,7 +180,7 @@ void BaseNRSolver<LinearSolver>::reset(){

template<class LinearSolver>
void BaseNRSolver<LinearSolver>::_dSbus_dV(const Eigen::Ref<const Eigen::SparseMatrix<cplx_type> > & Ybus,
const Eigen::Ref<const CplxVect > & V){
const Eigen::Ref<const CplxVect > & V){
auto timer = CustTimer();
const auto size_dS = V.size();
const CplxVect Vnorm = V.array() / V.array().abs();
Expand Down Expand Up @@ -331,6 +346,7 @@ void BaseNRSolver<LinearSolver>::fill_jacobian_matrix(const Eigen::SparseMatrix<
#endif // __COUT_TIMES
// first time i initialized the matrix, so i need to compute its sparsity pattern
fill_jacobian_matrix_unkown_sparsity_pattern(Ybus, V, slack_bus_id, slack_weights, pq, pvpq, pq_inv, pvpq_inv);
fill_value_map(slack_bus_id, pq, pvpq);
#ifdef __COUT_TIMES
std::cout << "\t\t fill_jacobian_matrix_unkown_sparsity_pattern : " << timer2.duration() << std::endl;
#endif // __COUT_TIMES
Expand All @@ -339,7 +355,8 @@ void BaseNRSolver<LinearSolver>::fill_jacobian_matrix(const Eigen::SparseMatrix<
// properly and faster (approx 3 times faster than the previous one)
#ifdef __COUT_TIMES
auto timer3 = CustTimer();
#endif // __COUT_TIMES
#endif //
if (BaseNRSolver<LinearSolver>::value_map_.size() == 0) fill_value_map(slack_bus_id,pq, pvpq);
fill_jacobian_matrix_kown_sparsity_pattern(slack_bus_id,
pq, pvpq
);
Expand Down Expand Up @@ -404,7 +421,7 @@ void BaseNRSolver<LinearSolver>::fill_jacobian_matrix_unkown_sparsity_pattern(
// optim : if the matrix was already computed, i don't initialize it, i instead reuse as much as i can
// i can do that because the matrix will ALWAYS have the same non zero coefficients.
// in this if, i allocate it in a "large enough" place to avoid copy when first filling it
if(J_.cols() != size_j) J_ = Eigen::SparseMatrix<real_type>(size_j,size_j);
if(J_.cols() != size_j) J_ = Eigen::SparseMatrix<real_type>(size_j, size_j);

std::vector<Eigen::Triplet<double> > coeffs; // HERE FOR PERF OPTIM (3)
coeffs.reserve(2*(dS_dVa_.nonZeros()+dS_dVm_.nonZeros()) + slack_weights.size()); // HERE FOR PERF OPTIM (3)
Expand Down Expand Up @@ -512,7 +529,6 @@ void BaseNRSolver<LinearSolver>::fill_jacobian_matrix_unkown_sparsity_pattern(
J_.setFromTriplets(coeffs.begin(), coeffs.end()); // HERE FOR PERF OPTIM (3)
// std::cout << "end fill jacobian unknown " << std::endl;
J_.makeCompressed();
fill_value_map(slack_bus_id, pq, pvpq);
// std::cout << "end fill_value_map" << std::endl;
}

Expand Down
75 changes: 50 additions & 25 deletions src/BaseNRSolverSingleSlack.tpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,15 +11,15 @@

template<class LinearSolver>
bool BaseNRSolverSingleSlack<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type> & Ybus,
CplxVect & V,
const CplxVect & Sbus,
const Eigen::VectorXi & slack_ids,
const RealVect & slack_weights, // unused here
const Eigen::VectorXi & pv,
const Eigen::VectorXi & pq,
int max_iter,
real_type tol
)
CplxVect & V,
const CplxVect & Sbus,
const Eigen::VectorXi & slack_ids,
const RealVect & slack_weights, // unused here
const Eigen::VectorXi & pv,
const Eigen::VectorXi & pq,
int max_iter,
real_type tol
)
{
/**
This method uses the newton raphson algorithm to compute voltage angles and magnitudes at each bus
Expand All @@ -31,18 +31,26 @@ bool BaseNRSolverSingleSlack<LinearSolver>::compute_pf(const Eigen::SparseMatrix
if(Sbus.size() != Ybus.rows() || Sbus.size() != Ybus.cols() ){
std::ostringstream exc_;
exc_ << "BaseNRSolverSingleSlack::compute_pf: Size of the Sbus should be the same as the size of Ybus. Currently: ";
exc_ << "Sbus (" << Sbus.size() << ") and Ybus (" << Ybus.rows() << ", " << Ybus.rows() << ").";
exc_ << "Sbus (" << Sbus.size() << ") and Ybus (" << Ybus.rows() << ", " << Ybus.cols() << ").";
throw std::runtime_error(exc_.str());
}
if(V.size() != Ybus.rows() || V.size() != Ybus.cols() ){
std::ostringstream exc_;
exc_ << "BaseNRSolverSingleSlack::compute_pf: Size of V (init voltages) should be the same as the size of Ybus. Currently: ";
exc_ << "V (" << V.size() << ") and Ybus (" << Ybus.rows()<<", "<<Ybus.rows() << ").";
exc_ << "V (" << V.size() << ") and Ybus (" << Ybus.rows()<<", "<<Ybus.cols() << ").";
throw std::runtime_error(exc_.str());
}
BaseNRSolver<LinearSolver>::reset_timer();

if(!BaseNRSolver<LinearSolver>::is_linear_solver_valid()) return false;
// std::cout << "singleslack" << std::endl;

if(BaseNRSolver<LinearSolver>::_solver_control.need_reset_solver() ||
BaseNRSolver<LinearSolver>::_solver_control.has_dimension_changed()){
BaseNRSolver<LinearSolver>::reset();
}

if(!BaseNRSolver<LinearSolver>::is_linear_solver_valid()){
return false;
}

BaseNRSolver<LinearSolver>::err_ = ErrorType::NoError; // reset the error if previous error happened
auto timer = CustTimer();
Expand Down Expand Up @@ -73,14 +81,20 @@ bool BaseNRSolverSingleSlack<LinearSolver>::compute_pf(const Eigen::SparseMatrix
bool has_just_been_initialized = false; // to avoid a call to klu_refactor follow a call to klu_factor in the same loop

const cplx_type m_i = BaseNRSolver<LinearSolver>::my_i; // otherwise it does not compile

BaseNRSolver<LinearSolver>::value_map_.clear(); // TODO smarter solver: only needed if ybus has changed or pq changed or pv changed
BaseNRSolver<LinearSolver>::dS_dVm_.resize(0,0); // TODO smarter solver: only needed if ybus has changed or pq changed or pv changed
BaseNRSolver<LinearSolver>::dS_dVa_.resize(0,0); // TODO smarter solver: only needed if ybus has changed or pq changed or pv changed
// BaseNRSolver<LinearSolver>::dS_dVm_.setZero(); // TODO smarter solver: only needed if ybus has changed
// BaseNRSolver<LinearSolver>::dS_dVa_.setZero(); // TODO smarter solver: only needed if ybus has changed
while ((!converged) & (BaseNRSolver<LinearSolver>::nr_iter_ < max_iter)){
BaseNRSolver<LinearSolver>::nr_iter_++;
// std::cout << "\tnr_iter_ " << BaseNRSolver<LinearSolver>::nr_iter_ << std::endl;
fill_jacobian_matrix(Ybus, BaseNRSolver<LinearSolver>::V_, pq, pvpq, pq_inv, pvpq_inv);
if(BaseNRSolver<LinearSolver>::need_factorize_){
BaseNRSolver<LinearSolver>::initialize();
if(BaseNRSolver<LinearSolver>::err_ != ErrorType::NoError){
// I got an error during the initialization of the linear system, i need to stop here
// std::cout << BaseNRSolver<LinearSolver>::err_ << std::endl;
res = false;
break;
}
Expand All @@ -95,6 +109,7 @@ bool BaseNRSolverSingleSlack<LinearSolver>::compute_pf(const Eigen::SparseMatrix
has_just_been_initialized = false;
if(BaseNRSolver<LinearSolver>::err_ != ErrorType::NoError){
// I got an error during the solving of the linear system, i need to stop here
// std::cout << BaseNRSolver<LinearSolver>::err_ << std::endl;
res = false;
break;
}
Expand All @@ -117,7 +132,11 @@ bool BaseNRSolverSingleSlack<LinearSolver>::compute_pf(const Eigen::SparseMatrix

F = BaseNRSolver<LinearSolver>::_evaluate_Fx(Ybus, BaseNRSolver<LinearSolver>::V_, Sbus, my_pv, pq);
bool tmp = F.allFinite();
if(!tmp) break; // divergence due to Nans
if(!tmp){
BaseNRSolver<LinearSolver>::err_ = ErrorType::InifiniteValue;
// std::cout << BaseNRSolver<LinearSolver>::err_ << std::endl;
break; // divergence due to Nans
}
converged = BaseNRSolver<LinearSolver>::_check_for_convergence(F, tol);
}
if(!converged){
Expand All @@ -135,17 +154,18 @@ bool BaseNRSolverSingleSlack<LinearSolver>::compute_pf(const Eigen::SparseMatrix
<< "\n\t timer_total_nr_: " << BaseNRSolver<LinearSolver>::timer_total_nr_
<< "\n\n";
#endif // __COUT_TIMES
BaseNRSolver<LinearSolver>::_solver_control.tell_none_changed();
return res;
}

template<class LinearSolver>
void BaseNRSolverSingleSlack<LinearSolver>::fill_jacobian_matrix(const Eigen::SparseMatrix<cplx_type> & Ybus,
const CplxVect & V,
const Eigen::VectorXi & pq,
const Eigen::VectorXi & pvpq,
const std::vector<int> & pq_inv,
const std::vector<int> & pvpq_inv
)
const CplxVect & V,
const Eigen::VectorXi & pq,
const Eigen::VectorXi & pvpq,
const std::vector<int> & pq_inv,
const std::vector<int> & pvpq_inv
)
{
/**
J has the shape
Expand All @@ -165,7 +185,6 @@ void BaseNRSolverSingleSlack<LinearSolver>::fill_jacobian_matrix(const Eigen::Sp
const int n_pvpq = static_cast<int>(pvpq.size());
const int n_pq = static_cast<int>(pq.size());
const int size_j = n_pvpq + n_pq;

// TODO to gain a bit more time below, try to compute directly, in _dSbus_dV(Ybus, V);
// TODO the `dS_dVa_[pvpq, pvpq]`
// TODO so that it's easier to retrieve in the next few lines !
Expand All @@ -177,6 +196,8 @@ void BaseNRSolverSingleSlack<LinearSolver>::fill_jacobian_matrix(const Eigen::Sp
#endif // __COUT_TIMES
// first time i initialized the matrix, so i need to compute its sparsity pattern
fill_jacobian_matrix_unkown_sparsity_pattern(Ybus, V, pq, pvpq, pq_inv, pvpq_inv);
fill_value_map(pq, pvpq);
// std::cout << "\t\tfill_jacobian_matrix_unkown_sparsity_pattern" << std::endl;
#ifdef __COUT_TIMES
std::cout << "\t\t fill_jacobian_matrix_unkown_sparsity_pattern : " << timer2.duration() << std::endl;
#endif // __COUT_TIMES
Expand All @@ -186,7 +207,12 @@ void BaseNRSolverSingleSlack<LinearSolver>::fill_jacobian_matrix(const Eigen::Sp
#ifdef __COUT_TIMES
auto timer3 = CustTimer();
#endif // __COUT_TIMES
if (BaseNRSolver<LinearSolver>::value_map_.size() == 0){
// std::cout << "\t\tfill_value_map called" << std::endl;
fill_value_map(pq, pvpq);
}
fill_jacobian_matrix_kown_sparsity_pattern(pq, pvpq);
// std::cout << "\t\tfill_jacobian_matrix_kown_sparsity_pattern" << std::endl;
#ifdef __COUT_TIMES
std::cout << "\t\t fill_jacobian_matrix_kown_sparsity_pattern : " << timer3.duration() << std::endl;
#endif // __COUT_TIMES
Expand Down Expand Up @@ -323,7 +349,6 @@ void BaseNRSolverSingleSlack<LinearSolver>::fill_jacobian_matrix_unkown_sparsity
}
// J_.setFromTriplets(coeffs.begin(), coeffs.end()); // HERE FOR PERF OPTIM (3)
BaseNRSolver<LinearSolver>::J_.makeCompressed();
fill_value_map(pq, pvpq);
}

/**
Expand All @@ -340,9 +365,9 @@ void BaseNRSolverSingleSlack<LinearSolver>::fill_value_map(
const int n_pvpq = static_cast<int>(pvpq.size());
BaseNRSolver<LinearSolver>::value_map_ = std::vector<cplx_type*> (BaseNRSolver<LinearSolver>::J_.nonZeros());

const int n_row = static_cast<int>(BaseNRSolver<LinearSolver>::J_.cols());
const int n_col = static_cast<int>(BaseNRSolver<LinearSolver>::J_.cols());
unsigned int pos_el = 0;
for (int col_=0; col_ < n_row; ++col_){
for (int col_=0; col_ < n_col; ++col_){
for (Eigen::SparseMatrix<real_type>::InnerIterator it(BaseNRSolver<LinearSolver>::J_, col_); it; ++it)
{
const int row_id = static_cast<int>(it.row());
Expand Down
1 change: 1 addition & 0 deletions src/BaseSolver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,6 +122,7 @@ bool BaseSolver::_check_for_convergence(const RealVect & F,
{
auto timer = CustTimer();
const auto norm_inf = F.lpNorm<Eigen::Infinity>();
// std::cout << "\t\tnorm_inf: " << norm_inf << std::endl;
bool res = norm_inf < tol;
timer_check_ += timer.duration();
return res;
Expand Down
11 changes: 10 additions & 1 deletion src/DCSolver.tpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,15 @@ bool BaseDCSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
// V is used the following way: at pq buses it's completely ignored. For pv bus only the magnitude is used,
// and for the slack bus both the magnitude and the angle are used.

if(!is_linear_solver_valid()) {
// err_ = ErrorType::NotInitError;
return false;
}
if(_solver_control.need_reset_solver() ||
_solver_control.has_dimension_changed()){
reset();
}

auto timer = CustTimer();
BaseSolver::reset_timer();
sizeYbus_with_slack_ = static_cast<int>(Ybus.rows());
Expand Down Expand Up @@ -135,7 +144,7 @@ bool BaseDCSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
#ifdef __COUT_TIMES
std::cout << "\t dc postproc: " << 1000. * timer_postproc.duration() << "ms" << std::endl;
#endif // __COUT_TIMES

_solver_control.tell_none_changed();
timer_total_nr_ += timer.duration();
return true;
}
Expand Down
Loading

0 comments on commit cc25e53

Please sign in to comment.