forked from Grid2op/lightsim2grid
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathBaseNRSolver.h
226 lines (199 loc) · 9.49 KB
/
BaseNRSolver.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
// Copyright (c) 2020, RTE (https://www.rte-france.com)
// See AUTHORS.txt
// This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
// If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
// you can obtain one at http://mozilla.org/MPL/2.0/.
// SPDX-License-Identifier: MPL-2.0
// This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
#ifndef BASENRSOLVER_H
#define BASENRSOLVER_H
#include "BaseSolver.h"
/**
Base class for Newton Raphson based solver
**/
template<class LinearSolver>
class BaseNRSolver : public BaseSolver
{
public:
BaseNRSolver():need_factorize_(true), timer_initialize_(0.), timer_dSbus_(0.), timer_fillJ_(0.) {}
virtual
Eigen::Ref<const Eigen::SparseMatrix<real_type> > get_J() const {
return J_;
}
virtual
Eigen::SparseMatrix<real_type> get_J_python() const {
Eigen::SparseMatrix<real_type> res = get_J();
return res;
}
virtual
std::tuple<double, double, double, double, double, double, double> get_timers_jacobian()
{
// TODO refacto that, and change the order
auto res = std::tuple<double, double, double, double, double, double, double>(
timer_Fx_, timer_solve_, timer_initialize_, timer_check_, timer_dSbus_, timer_fillJ_, timer_total_nr_);
return res;
}
virtual
bool 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
) ;
virtual void reset();
protected:
virtual void reset_timer(){
BaseSolver::reset_timer();
timer_dSbus_ = 0.;
timer_fillJ_ = 0.;
timer_initialize_ = 0.;
}
virtual
void initialize(){
auto timer = CustTimer();
n_ = static_cast<int>(J_.cols()); // should be equal to J_.nrows()
err_ = ErrorType::NoError; // reset error message
const ErrorType init_status = _linear_solver.initialize(J_);
if(init_status != ErrorType::NoError){
// std::cout << "init_ok " << init_ok << std::endl;
err_ = init_status;
}
need_factorize_ = false;
timer_initialize_ += timer.duration();
}
virtual
void solve(RealVect & b, bool has_just_been_inialized){
auto timer = CustTimer();
const ErrorType solve_status = _linear_solver.solve(J_, b, has_just_been_inialized);
if(solve_status != ErrorType::NoError){
// std::cout << "solve error: " << solve_status << std::endl;
err_ = solve_status;
}
timer_solve_ += timer.duration();
}
void _dSbus_dV(const Eigen::Ref<const Eigen::SparseMatrix<cplx_type> > & Ybus,
const Eigen::Ref<const CplxVect > & V);
void _get_values_J(int & nb_obj_this_col,
std::vector<Eigen::Index> & inner_index,
std::vector<real_type> & values,
const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & mat, // ex. dS_dVa_r
const std::vector<int> & index_row_inv, // ex. pvpq_inv
const Eigen::VectorXi & index_col, // ex. pvpq
Eigen::Index col_id,
Eigen::Index row_lag, // 0 for J11 for example, n_pvpq for J12
Eigen::Index col_lag
);
void _get_values_J(int & nb_obj_this_col,
std::vector<Eigen::Index> & inner_index,
std::vector<real_type> & values,
const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & mat, // ex. dS_dVa_r
const std::vector<int> & index_row_inv, // ex. pvpq_inv
Eigen::Index col_id_mat, // ex. pvpq(col_id)
Eigen::Index row_lag, // 0 for J11 for example, n_pvpq for J12
Eigen::Index col_lag // to remove the ref slack bus from this
);
void fill_jacobian_matrix(const Eigen::SparseMatrix<cplx_type> & Ybus,
const CplxVect & V,
Eigen::Index slack_bus_id,
const RealVect & slack_weights,
const Eigen::VectorXi & pq,
const Eigen::VectorXi & pvpq,
const std::vector<int> & pq_inv,
const std::vector<int> & pvpq_inv
);
void fill_jacobian_matrix_kown_sparsity_pattern(
Eigen::Index slack_bus_id,
const Eigen::VectorXi & pq,
const Eigen::VectorXi & pvpq
);
void fill_jacobian_matrix_unkown_sparsity_pattern(
const Eigen::SparseMatrix<cplx_type> & Ybus,
const CplxVect & V,
Eigen::Index slack_bus_id,
const RealVect & slack_weights,
const Eigen::VectorXi & pq,
const Eigen::VectorXi & pvpq,
const std::vector<int> & pq_inv,
const std::vector<int> & pvpq_inv
);
void fill_value_map(Eigen::Index slack_bus_id,
const Eigen::VectorXi & pq,
const Eigen::VectorXi & pvpq);
protected:
// used linear solver
LinearSolver _linear_solver;
// solution of the problem
Eigen::SparseMatrix<real_type> J_; // the jacobian matrix
Eigen::SparseMatrix<cplx_type> dS_dVm_;
Eigen::SparseMatrix<cplx_type> dS_dVa_;
bool need_factorize_;
// to store the mapping from the element of J_ in dS_dVm_ and dS_dVa_
// it does not own any memory at all !
std::vector<cplx_type*> value_map_;
// timers
double timer_initialize_;
double timer_dSbus_;
double timer_fillJ_;
Eigen::SparseMatrix<real_type>
create_jacobian_matrix_test(const Eigen::SparseMatrix<cplx_type> & Ybus,
const CplxVect & V,
const RealVect & slack_weights,
const Eigen::VectorXi & pq,
const Eigen::VectorXi & pvpq
){
// DO NOT USE, FOR DEBUG ONLY (especially for multiple slacks)
const auto & n_pvpq = pvpq.size();
const auto & n_pq = pvpq.size();
std::vector<int> pvpq_inv(V.size(), -1);
for(int inv_id=0; inv_id < n_pvpq; ++inv_id) pvpq_inv[pvpq(inv_id)] = inv_id;
std::vector<int> pq_inv(V.size(), -1);
for(int inv_id=0; inv_id < n_pq; ++inv_id) pq_inv[pq(inv_id)] = inv_id;
// TODO if bug when using it, check the "pvpq" below,
// in theory its "pv" !
int slack_bus_id = extract_slack_bus_id(pvpq, pq,
static_cast<unsigned int>(V.size())
);
fill_jacobian_matrix(Ybus, V, static_cast<Eigen::Index>(slack_bus_id),
slack_weights, pq, pvpq, pq_inv, pvpq_inv);
return J_;
}
private:
// no copy allowed
BaseNRSolver( const BaseNRSolver & ) =delete ;
BaseNRSolver & operator=( const BaseNRSolver & ) =delete ;
/** helper function to print the max_col left most columns of the J matrix **/
void print_J(int min_col=-1, int max_col=-1) const{
auto size_J = J_.cols();
if(max_col == -1) max_col = static_cast<int>(size_J);
if(min_col == -1) min_col = 0;
for (int col_id=min_col; col_id < max_col; ++col_id){
for (Eigen::SparseMatrix<real_type>::InnerIterator it(J_, col_id); it; ++it)
{
std::cout << it.row() << ", " << it.col() << ": " << it.value() << std::endl;
}
std::cout << std::endl;
}
}
private:
Eigen::SparseMatrix<cplx_type>
_make_diagonal_matrix(const Eigen::Ref<const CplxVect > & diag_val){
// TODO their might be a more efficient way to do that
auto n = diag_val.size();
Eigen::SparseMatrix<cplx_type> res(n,n);
// first method, without a loop of mine
res.setIdentity(); // segfault if attempt to use this function without this
res.diagonal() = diag_val;
// second method, with an "optimized" loop
// res.reserve(Eigen::VectorXi::Constant(n,1)); // i reserve one number per columns (speed optim)
// for(unsigned int i = 0; i < n; ++i){
// res.insert(i,i) = diag_val(i);
//}
return res;
}
};
#include "BaseNRSolver.tpp"
#endif // BASENRSOLVER_H