Skip to content

Commit 8b003c8

Browse files
mgoversmgovers
authored
Merge pull request #1150 from PowerGridModel/feature/statistics-2nd-order-variance-approximation
Current sensor statistics: to polar to decomposed statistics
2 parents c47046b + 29e11ed commit 8b003c8

File tree

3 files changed

+206
-78
lines changed

3 files changed

+206
-78
lines changed

‎power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp‎

Lines changed: 71 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -126,6 +126,41 @@ template <symmetry_tag sym_type> struct DecomposedComplexRandVar {
126126
}
127127
};
128128

129+
namespace detail {
130+
// Var(f_i) ≈ \Sum_j Var(x_j) * ||df_i/dx_j||^2 // first order
131+
// + 1/2 * \Sum_{jk} Var(x_j) * ||d^2 f_i/dx_jdx_k||^2 Var(x_k) // second order
132+
//
133+
// Re{I} = I cos(θ)
134+
// Im{I} = I sin(θ)
135+
//
136+
// dRe{I}/dI = cos(θ) dIm{I}/dI = sin(θ)
137+
// dRe{I}/dθ = -I sin(θ) dIm{I}/dθ = I cos(θ)
138+
// d^2 Re{I}/dI^2 = 0 d^2 Im{I}/dI^2 = 0
139+
// d^2 Re{I}/dθ^2 = -I cos(θ) d^2 Im{I}/dθ^2 = -I sin(θ)
140+
// d^2 Re{I}/dIdθ = -sin(θ) d^2 Im{I}/dIdθ = cos(θ)
141+
//
142+
// Var(Re{I}) = Var(I cos(θ))
143+
// ≈ Var(I) * cos^2(θ) + Var(θ) * I^2 * sin^2(θ) // first order
144+
// + 1/2 * Var(θ)^2 * I^2 * cos^2(θ) + Var(I) * Var(θ) * sin^2(θ) // second order
145+
// Var(Im{I}) = Var(I sin(θ))
146+
// ≈ Var(I) * sin^2(θ) + Var(θ) * I^2 * cos^2(θ) // first order
147+
// + 1/2 * Var(θ)^2 * I^2 * sin^2(θ) + Var(I) * Var(θ) * cos^2(θ) // second order
148+
inline auto compute_decomposed_variance_from_polar(auto const& magnitude, auto const& cos_angle, auto const& sin_angle,
149+
auto const& magnitude_variance, auto const& angle_variance) {
150+
auto const cos2_angle = cos_angle * cos_angle;
151+
auto const sin2_angle = sin_angle * sin_angle;
152+
auto const magnitude2 = magnitude * magnitude;
153+
154+
// second order approximation
155+
return std::make_pair(magnitude_variance * cos2_angle + magnitude2 * angle_variance * sin2_angle +
156+
0.5 * magnitude2 * angle_variance * angle_variance * cos2_angle +
157+
magnitude_variance * angle_variance * sin2_angle,
158+
magnitude_variance * sin2_angle + magnitude2 * angle_variance * cos2_angle +
159+
0.5 * magnitude2 * angle_variance * angle_variance * sin2_angle +
160+
magnitude_variance * angle_variance * cos2_angle);
161+
}
162+
} // namespace detail
163+
129164
// Complex measured value of a sensor in polar coordinates (magnitude and angle)
130165
// (rotationally symmetric)
131166
template <symmetry_tag sym_type> struct PolarComplexRandVar {
@@ -145,65 +180,66 @@ template <symmetry_tag sym_type> struct PolarComplexRandVar {
145180
}
146181

147182
// For sym to sym conversion:
148-
// Var(I_Re) ≈ Var(I) * cos^2(θ) + Var(θ) * I^2 * sin^2(θ)
149-
// Var(I_Im) ≈ Var(I) * sin^2(θ) + Var(θ) * I^2 * cos^2(θ)
183+
// Var(I_Re) ≈ Var(I) * cos^2(θ) + Var(θ) * I^2 * sin^2(θ) // first order
184+
// + 1/2 * Var(θ)^2 * I^2 * cos^2(θ) + Var(I) * Var(θ) * sin^2(θ) // second order
185+
// Var(I_Im) ≈ Var(I) * sin^2(θ) + Var(θ) * I^2 * cos^2(θ) // first order
186+
// + 1/2 * Var(θ)^2 * I^2 * sin^2(θ) + Var(I) * Var(θ) * cos^2(θ) // second order
150187
// For asym to asym conversion:
151-
// Var(I_Re,p) ≈ Var(I_p) * cos^2(θ_p) + Var(θ_p) * I_p^2 * sin^2(θ_p)
152-
// Var(I_Im,p) ≈ Var(I_p) * sin^2(θ_p) + Var(θ_p) * I_p^2 * cos^2(θ_p)
188+
// Var(I_Re,p) ≈ Var(I_p) * cos^2(θ_p) + Var(θ_p) * I_p^2 * sin^2(θ_p) + 2nd order terms (idem)
189+
// Var(I_Im,p) ≈ Var(I_p) * sin^2(θ_p) + Var(θ_p) * I_p^2 * cos^2(θ_p) + 2nd order terms (idem)
153190
explicit operator DecomposedComplexRandVar<sym>() const {
154191
auto const cos_theta = cos(angle.value);
155192
auto const sin_theta = sin(angle.value);
156193
auto const real_component = magnitude.value * cos_theta;
157194
auto const imag_component = magnitude.value * sin_theta;
158-
return DecomposedComplexRandVar<sym>{
159-
.real_component = {.value = real_component,
160-
.variance = magnitude.variance * cos_theta * cos_theta +
161-
imag_component * imag_component * angle.variance},
162-
.imag_component = {.value = imag_component,
163-
.variance = magnitude.variance * sin_theta * sin_theta +
164-
real_component * real_component * angle.variance}};
195+
196+
auto const [real_variance, imag_variance] = detail::compute_decomposed_variance_from_polar(
197+
magnitude.value, cos_theta, sin_theta, magnitude.variance, angle.variance);
198+
199+
return DecomposedComplexRandVar<sym>{.real_component = {.value = real_component, .variance = real_variance},
200+
.imag_component = {.value = imag_component, .variance = imag_variance}};
165201
}
166202

167-
// Var(I_Re,p) ≈ Var(I) * cos^2(θ - 2πp/3) + Var(θ) * I^2 * sin^2(θ - 2πp/3)
168-
// Var(I_Im,p) ≈ Var(I) * sin^2(θ - 2πp/3) + Var(θ) * I^2 * cos^2(θ - 2πp/3)
203+
// Var(I_Re,p) ≈ Var(I) * cos^2(θ - 2πp/3) + Var(θ) * I^2 * sin^2(θ - 2πp/3) + 2nd order terms
204+
// Var(I_Im,p) ≈ Var(I) * sin^2(θ - 2πp/3) + Var(θ) * I^2 * cos^2(θ - 2πp/3) + 2nd order terms
169205
explicit operator DecomposedComplexRandVar<asymmetric_t>() const
170206
requires(is_symmetric_v<sym>)
171207
{
172-
ComplexValue<asymmetric_t> const unit_complex{exp(1.0i * angle.value)};
208+
ComplexValue<asymmetric_t> const unit_complex{exp(1.0i * angle.value)}; // correctly handle phase shifts
173209
ComplexValue<asymmetric_t> const complex = magnitude.value * unit_complex;
210+
RealValue<asymmetric_t> const cos_theta = real(unit_complex);
211+
RealValue<asymmetric_t> const sin_theta = imag(unit_complex);
212+
RealValue<asymmetric_t> const real_component = real(complex);
213+
RealValue<asymmetric_t> const imag_component = imag(complex);
214+
215+
auto const [real_variance, imag_variance] = detail::compute_decomposed_variance_from_polar(
216+
magnitude.value, cos_theta, sin_theta, magnitude.variance, angle.variance);
217+
174218
return DecomposedComplexRandVar<asymmetric_t>{
175-
.real_component = {.value = real(complex),
176-
.variance = magnitude.variance * real(unit_complex) * real(unit_complex) +
177-
imag(complex) * imag(complex) * angle.variance},
178-
.imag_component = {.value = imag(complex),
179-
.variance = magnitude.variance * imag(unit_complex) * imag(unit_complex) +
180-
real(complex) * real(complex) * angle.variance}};
219+
.real_component = {.value = real_component, .variance = real_variance},
220+
.imag_component = {.value = imag_component, .variance = imag_variance}};
181221
}
182222

183223
// Var(I_Re) ≈ (1 / 9) * sum_p(Var(I_p) * cos^2(theta_p + 2 * pi * p / 3) + Var(theta_p) * I_p^2 * sin^2(theta_p + 2
184-
// * pi * p / 3))
224+
// * pi * p / 3)) + 2nd order terms
185225
// Var(I_Im) ≈ (1 / 9) * sum_p(Var(I_p) * sin^2(theta_p + 2 * pi * p / 3) + Var(theta_p) * I_p^2 *
186-
// cos^2(theta_p + 2 * pi * p / 3))
226+
// cos^2(theta_p + 2 * pi * p / 3)) + 2nd order terms
187227
explicit operator DecomposedComplexRandVar<symmetric_t>() const
188228
requires(is_asymmetric_v<sym>)
189229
{
190230
ComplexValue<asymmetric_t> const unit_complex{exp(1.0i * angle.value)};
191231
ComplexValue<asymmetric_t> const unit_pos_seq_per_phase{unit_complex(0), a * unit_complex(1),
192232
a2 * unit_complex(2)};
193233
DoubleComplex const pos_seq_value = pos_seq(magnitude.value * unit_complex);
234+
RealValue<asymmetric_t> const cos_theta = real(unit_pos_seq_per_phase);
235+
RealValue<asymmetric_t> const sin_theta = imag(unit_pos_seq_per_phase);
236+
237+
auto const [real_variance, imag_variance] = detail::compute_decomposed_variance_from_polar(
238+
magnitude.value, cos_theta, sin_theta, magnitude.variance, angle.variance);
239+
194240
return DecomposedComplexRandVar<symmetric_t>{
195-
.real_component = {.value = real(pos_seq_value),
196-
.variance = sum_val(magnitude.variance * real(unit_pos_seq_per_phase) *
197-
real(unit_pos_seq_per_phase) +
198-
imag(unit_pos_seq_per_phase) * imag(unit_pos_seq_per_phase) *
199-
magnitude.value * magnitude.value * angle.variance) /
200-
9.0},
201-
.imag_component = {
202-
.value = imag(pos_seq_value),
203-
.variance = sum_val(magnitude.variance * imag(unit_pos_seq_per_phase) * imag(unit_pos_seq_per_phase) +
204-
real(unit_pos_seq_per_phase) * real(unit_pos_seq_per_phase) * magnitude.value *
205-
magnitude.value * angle.variance) /
206-
9.0}};
241+
.real_component = {.value = real(pos_seq_value), .variance = sum_val(real_variance) / 9.0},
242+
.imag_component = {.value = imag(pos_seq_value), .variance = sum_val(imag_variance) / 9.0}};
207243
}
208244
};
209245

@@ -226,7 +262,7 @@ template <symmetry_tag sym> inline auto conj(PolarComplexRandVar<sym> var) {
226262
}
227263

228264
namespace statistics {
229-
// Var(s x) ≈ Var(x) * ||s||²
265+
// Var(s x) ≈ Var(x) * ||s||^2
230266
template <symmetry_tag sym, template <symmetry_tag> typename RandVarType>
231267
requires is_in_list_c<RandVarType<sym>, UniformRealRandVar<sym>, IndependentRealRandVar<sym>,
232268
UniformComplexRandVar<sym>, IndependentComplexRandVar<sym>>

‎tests/cpp_unit_tests/test_current_sensor.cpp‎

Lines changed: 26 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -82,11 +82,17 @@ TEST_CASE("Test current sensor") {
8282
// Check symmetric sensor output for symmetric parameters
8383
CHECK(sym_sensor_param.angle_measurement_type == angle_measurement_type);
8484
// Var(I_Re) ≈ Var(I) * cos^2(pi/4) + Var(θ) * I^2 * sin^2(pi/4)
85+
// + (1/2) * Var(θ)^2 * I^2 * cos^2(pi/4) + Var(I) * Var(θ) * sin^2(pi/4)
8586
CHECK(sym_sensor_param.measurement.real_component.variance ==
86-
doctest::Approx(0.5 * (i_variance_pu + i_angle_variance_pu * i_pu * i_pu)));
87+
doctest::Approx(0.5 * (i_variance_pu + i_angle_variance_pu * i_pu * i_pu +
88+
0.5 * i_angle_variance_pu * i_angle_variance_pu * i_pu * i_pu +
89+
i_variance_pu * i_angle_variance_pu)));
8790
// Var(I_Im) ≈ Var(I) * sin^2(pi/4) + Var(θ) * I^2 * cos^2(pi/4)
91+
// + (1/2) * Var(θ)^2 * I^2 * sin^2(pi/4) + Var(I) * Var(θ) * cos^2(pi/4)
8892
CHECK(sym_sensor_param.measurement.imag_component.variance ==
89-
doctest::Approx(0.5 * (i_variance_pu + i_angle_variance_pu * i_pu * i_pu)));
93+
doctest::Approx(0.5 * (i_variance_pu + i_angle_variance_pu * i_pu * i_pu +
94+
0.5 * i_angle_variance_pu * i_angle_variance_pu * i_pu * i_pu +
95+
i_variance_pu * i_angle_variance_pu)));
9096
CHECK(real(sym_sensor_param.measurement.value()) == doctest::Approx(i_pu * cos(i_angle)));
9197
CHECK(imag(sym_sensor_param.measurement.value()) == doctest::Approx(i_pu * sin(i_angle)));
9298

@@ -114,11 +120,17 @@ TEST_CASE("Test current sensor") {
114120
CHECK(asym_sensor_param.angle_measurement_type == angle_measurement_type);
115121

116122
CHECK(asym_sensor_param.measurement.real_component.variance[0] ==
117-
doctest::Approx(0.5 * (i_variance_pu + i_angle_variance_pu * i_pu * i_pu)));
123+
doctest::Approx(0.5 * (i_variance_pu + i_angle_variance_pu * i_pu * i_pu +
124+
0.5 * i_angle_variance_pu * i_angle_variance_pu * i_pu * i_pu +
125+
i_variance_pu * i_angle_variance_pu)));
118126
auto const shifted_i_angle = i_angle + deg_240;
119-
CHECK(asym_sensor_param.measurement.imag_component.variance[1] ==
120-
doctest::Approx(i_variance_pu * sin(shifted_i_angle) * sin(shifted_i_angle) +
121-
i_angle_variance_pu * i_pu * i_pu * cos(shifted_i_angle) * cos(shifted_i_angle)));
127+
CHECK(
128+
asym_sensor_param.measurement.imag_component.variance[1] ==
129+
doctest::Approx(i_variance_pu * sin(shifted_i_angle) * sin(shifted_i_angle) +
130+
i_angle_variance_pu * i_pu * i_pu * cos(shifted_i_angle) * cos(shifted_i_angle) +
131+
0.5 * i_angle_variance_pu * i_angle_variance_pu * i_pu * i_pu *
132+
sin(shifted_i_angle) * sin(shifted_i_angle) +
133+
i_variance_pu * i_angle_variance_pu * cos(shifted_i_angle) * cos(shifted_i_angle)));
122134
CHECK(real(asym_sensor_param.measurement.value()[0]) == doctest::Approx(i_pu * cos(i_angle)));
123135
CHECK(imag(asym_sensor_param.measurement.value()[1]) == doctest::Approx(i_pu * sin(shifted_i_angle)));
124136

@@ -162,8 +174,10 @@ TEST_CASE("Test current sensor") {
162174
auto const sym_param = sym_current_sensor.calc_param<symmetric_t>();
163175

164176
CHECK(sym_param.angle_measurement_type == angle_measurement_type);
165-
CHECK(sym_param.measurement.real_component.variance == doctest::Approx(pow(1.0 / base_current, 2)));
166-
CHECK(sym_param.measurement.imag_component.variance == doctest::Approx(pow(0.2 / base_current, 2)));
177+
CHECK(sym_param.measurement.real_component.variance ==
178+
doctest::Approx(pow(1.0 / base_current, 2) * (1.0 + 0.5 * pow(0.2, 4))));
179+
CHECK(sym_param.measurement.imag_component.variance ==
180+
doctest::Approx(pow(1.0 / base_current, 2) * (pow(0.2, 2) + pow(1.0 * 0.2, 2))));
167181
CHECK(real(sym_param.measurement.value()) == doctest::Approx(1.0 / base_current));
168182
CHECK(imag(sym_param.measurement.value()) == doctest::Approx(0.0 / base_current));
169183
}
@@ -174,8 +188,10 @@ TEST_CASE("Test current sensor") {
174188
auto const sym_param = sym_current_sensor.calc_param<symmetric_t>();
175189

176190
CHECK(sym_param.angle_measurement_type == angle_measurement_type);
177-
CHECK(sym_param.measurement.real_component.variance == doctest::Approx(pow(0.2 / base_current, 2)));
178-
CHECK(sym_param.measurement.imag_component.variance == doctest::Approx(pow(1.0 / base_current, 2)));
191+
CHECK(sym_param.measurement.real_component.variance ==
192+
doctest::Approx(pow(1.0 / base_current, 2) * (pow(0.2, 2) + pow(0.2, 2))));
193+
CHECK(sym_param.measurement.imag_component.variance ==
194+
doctest::Approx(pow(1.0 / base_current, 2) * (1.0 + 0.5 * pow(0.2, 4))));
179195
CHECK(real(sym_param.measurement.value()) == doctest::Approx(0.0 / base_current));
180196
CHECK(imag(sym_param.measurement.value()) == doctest::Approx(1.0 / base_current));
181197
}

0 commit comments

Comments
 (0)