[FEATURE] Support current measurement in state estimation

[TonyXiang8787]

Background

Currently the state estimation supports two types of sensors: voltage_sensor and power_sensor. Their meaning is self explanatory. In distribution grids, we have also many current measurements. Supporting current magnitude measurement is difficult since we do not know the direction of active/reactive power. To decouple the current, we need the voltage as a reference, with which we can calculate power. Therefore, we only support power measurements as this moment.

Current measurement with relative phase angle to the local voltage

Sometimes we have both voltage and current measurements in practice. However, certain voltage magnitude measurement has prohibitable high error margin which makes the calculated power value useless. Meanwhile, the phase angle between the voltage and current is actually accurately measured. In this case, we have current magnitude + relative phase angle measurement.

Current measurement using PMU with global phase angle

If we install PMU current measurement in the grid, we can measure the current magnitude and the global phase angle with regarding the reference angle in the HV/MV substation. In this way we measure the current phasor.

Feature request

This issue proposes to support decoupled current measurement into PGM state estimation. Concretely, the following tasks:

• Define new component sym_current_sensor and asym_current_sensor.
  • They should have the attributes i_measured and i_angle_measured, representing the current magnitude and angle.
  • They should have the attributes angle_measurement_type which is an enumeration of local and global.
    • In case of local, this is a relative angle to the local voltage.
    • In case of global, this is a PMU current phasor with the global angle to the reference angle in the HV/MV substation.
  • The reference direction should be the same as the measured object, similar to power_sensor.
  • Also define relevant sigmas.
  • They can only measure branch and branch3
• In the observability check, treat the current sensor the same as power sensor. They both contribute to observability.
• In iterative linear solver, the implementation is straightforward, as the measured value in linear solver is already current phasor. In case of a local angle, we just need to shift the current phase angle with the node voltage phase angle.
  • Mathematics are refined, reviewed and approved in the below section
• In Newton-Raphson solver, we need to refine the mathematics of current sensor and implement them.
  • Mathematics are refined, reviewed and approved in the below section

Mathematics

In all variance equations below, the following approximation is used that holds for both real and complex functions, as well as complex and vector-valued functions of multiple parameters

$$\begin{eqnarray} \text{let } \vec{F}(\vec{x}) \equiv \sum_i \vec{e}_i f_i(x_j) \\\ \text{then } \text{Var}(f_i) \approx \sum_j \text{Var}(x_j) \|\frac{\delta f_i}{\delta{x_j}}(\vec{x})\|^2 \end{eqnarray}$$

This allows us to treat the imaginary and real components separately, and then combine the two again using $$\text{Var}(z) = \text{Var}(\Re(z)) + \text{Var}(\Im(z))$$, which in turn also is a special case of the above equation.

Symmetric current sensors

As intermediate calculation parameters, we model the independent measurements of the magnitude and angle to a single complex measurement with (assumed to be) independent real and imaginary components. Note that this is an approximation. The variances are calculated as follows. (cfr. #861 (comment) )

Symmetric calculations

$$\begin{eqnarray} && && \Re\{I e^{i\theta}\} &&=&& I \cos \theta \\\ && && \Im\{I e^{i\theta}\} &&=&& I \sin \theta \\\ \text{Var}(\Re\{I e^{i\theta}\}) &&=&& \text{Var}(I \cos \theta) &&\approx && \text{Var}(I)\cos^2 \theta + I^2 \text{Var}(\theta) \sin^2 \theta \\\ \text{Var}(\Im\{I e^{i\theta}\}) &&=&& \text{Var}(I \sin \theta) &&\approx && \text{Var}(I)\sin^2 \theta + I^2 \text{Var}(\theta) \cos^2 \theta \\\ \text{Var}(I e^{i\theta}) &&=&& \text{Var}(\Re\{I e^{i\theta}\}) + \text{Var}(\Im\{I e^{i\theta}\}) && \approx && \text{Var}(I) + I^2 \text{Var}(\theta) \\\ \sigma_{\Re} &&=&& \sqrt{\text{Var}(\Re\{I e^{i\theta}\})} && \approx && \sqrt{\text{Var}(I)\cos^2 \theta + I^2 \text{Var}(\theta) \sin^2 \theta} \\\ \sigma_{\Im} &&=&& \sqrt{\text{Var}(\Im\{I e^{i\theta}\})} && \approx && \sqrt{\text{Var}(I)\sin^2 \theta + I^2 \text{Var}(\theta) \cos^2 \theta} \\\ \sigma_{\underline{I}} &&=&& \sqrt{\text{Var}(I e^{i\theta})} \\\ &&=&& \sqrt{\sigma_{\Im}^2 + \sigma_{\Re}^2} && \approx && \sqrt{\text{Var}(I) + I^2 \text{Var}(\theta)} \end{eqnarray}$$

This amounts to a box that weighs the extrapolated error bars to obtain the corners.

• For NRSE, we can use $$\sigma_{\Re}$$ and $$\sigma_{\Im}$$ directly.
• For ILSE, we can use $$\sigma_{\underline{I}} = \sqrt{\sigma_{\Im}^2 + \sigma_{\Re}^2}$$

Asymmetric calculations

Let $$p\in{a,b,c}$$ be the phase index

Global angle measurement

For global angle measurements, the following needs to be done:

$$\begin{eqnarray} a &&=&& e^{2 \pi i/3} \\\ \underline{I}_p &&=&& a^{-p} I e^{i \theta} && = && I e^{i \theta - \frac{2\pi i}{3} p} \\\ I_{\Re,p} &&\coloneqq && \Re\{\underline{I}_p\} &&=&& I \cos (\theta - \frac{2\pi}{3}p) \\\ I_{\Im,p} &&\coloneqq && \Im\{\underline{I}_p\} &&=&& I \sin (\theta - \frac{2\pi}{3}p) \\\ && && \text{Var}(I_{\Re,p}) &&\approx && \text{Var}(I) \cos^2 (\theta - \frac{2\pi}{3}p) + \text{Var}(\theta)I^2 \sin^2 (\theta - \frac{2\pi}{3}p) \\\ && && \text{Var}(I_{\Im,p}) &&\approx && \text{Var}(I) \sin^2 (\theta - \frac{2\pi}{3}p) + \text{Var}(\theta)I^2 \cos^2 (\theta - \frac{2\pi}{3}p) \\\ && && \text{Var}(\underline{I}_p) &&\approx && \text{Var}(I) + \text{Var}(\theta)I^2 \\\ \sigma_{\Re,p} &&=&& \sqrt{\text{Var}(\Re\{I e^{i\theta - \frac{2\pi i}{3}p}\})} && \approx && \sqrt{\text{Var}(I)\cos^2 (\theta - \frac{2\pi}{3}p) + I^2 \text{Var}(\theta) \sin^2 (\theta - \frac{2\pi}{3}p)} \\\ \sigma_{\Im,p} &&=&& \sqrt{\text{Var}(\Im\{I e^{i\theta - \frac{2\pi i}{3}p}\})} && \approx && \sqrt{\text{Var}(I)\sin^2 (\theta - \frac{2\pi}{3}p) + I^2 \text{Var}(\theta) \cos^2 (\theta - \frac{2\pi}{3}p)} \\\ \sigma_{p} &&=&& \sqrt{\text{Var}(I e^{i\theta - \frac{2\pi i}{3}p})} \\\ &&=&& \sqrt{\sigma_{\Im,p}^2 + \sigma_{\Re,p}^2} && \approx && \sqrt{\text{Var}(I) + I^2 \text{Var}(\theta)} \\\ && && && \equiv && \sigma \end{eqnarray}$$

Local angle measurement

For local angle measurement, no phase shift should be introduced

$$\begin{eqnarray} \underline{I}_p &&=&& I e^{i \theta} && = && I e^{i \theta} \\\ I_{\Re,p} &&\coloneqq && \Re\{\underline{I}_p\} &&=&& I \cos \theta \\\ I_{\Im,p} &&\coloneqq && \Im\{\underline{I}_p\} &&=&& I \sin \theta \\\ && && \text{Var}(I_{\Re,p}) &&\approx && \text{Var}(I) \cos^2 \theta + \text{Var}(\theta)I^2 \sin^2 \theta \\\ && && \text{Var}(I_{\Im,p}) &&\approx && \text{Var}(I) \sin^2 \theta + \text{Var}(\theta)I^2 \cos^2 \theta \\\ && && \text{Var}(\underline{I}_p) &&\approx && \text{Var}(I) + \text{Var}(\theta)I^2 \\\ \sigma_{\Re,p} &&=&& \sqrt{\text{Var}(\Re\{I e^{i\theta}\})} && \approx && \sqrt{\text{Var}(I)\cos^2 \theta + I^2 \text{Var}(\theta) \sin^2 \theta} \\\ && && && \equiv && \sigma_{\Re} \\\ \sigma_{\Im,p} &&=&& \sqrt{\text{Var}(\Im\{I e^{i\theta}\})} && \approx && \sqrt{\text{Var}(I)\sin^2 \theta + I^2 \text{Var}(\theta) \cos^2 \theta} \\\ && && && \equiv && \sigma_{\Im} \\\ \sigma_{p} &&=&& \sqrt{\text{Var}(I e^{i\theta})} \\\ &&=&& \sqrt{\sigma_{\Im,p}^2 + \sigma_{\Re,p}^2} && \approx && \sqrt{\text{Var}(I) + I^2 \text{Var}(\theta)} \\\ && && && \equiv && \sigma \end{eqnarray}$$

Asymmetric calculations on symmetric sensors with local current measurements are, therefore, component-wise the same as for symmetric calculations.

Asymmetric current sensors

The measurements are $$I_p$$ and $$\theta_p$$ for phases $$p$$
We assume the variance for the magnitude and the angle to be the same for all phases: $$\text{Var}(I_p) \equiv \text{Var}(I_m)$$ and $$\text{Var}(\theta_p)\equiv \text{Var}(\theta_m)$$, $$\forall p$$

Symmetric calculations

Global angle measurement

For global angle measurements, the symmetric current and phase are those of the positive sequence current.

$$\begin{eqnarray} a &&=&& e^{2 \pi i/3} \\\ \underline{I} \coloneqq \underline{I}_1 &&=&& \frac{\sum_p a^p I_p e^{i\theta_p}}{3} \\\ I &&=&& \|\underline{I}\| \\\ e^{i \theta} &&=&& \frac{\underline{I}}{I} \\\ \theta &&=&& \arctan(\frac{\Im\{e^{i\theta}\}}{\Re\{e^{i\theta}\}}) \\\ \text{Var}(I_{\Re}) &&\approx && \frac{1}{9} \sum_p(\text{Var}(I_p) \cos^2(\theta_p + \frac{2 \pi}{3}p) + \text{Var}(\theta_p)I_p^2 \sin^2(\theta_p + \frac{2\pi}{3}p)) \\\ \text{Var}(I_{\Im}) &&\approx && \frac{1}{9} \sum_p(\text{Var}(I_p) \sin^2(\theta_p + \frac{2 \pi}{3}p) + \text{Var}(\theta_p)I_p^2 \cos^2(\theta_p + \frac{2\pi}{3}p)) \\\ \text{Var}(\underline{I}) &&\approx && \frac{1}{9} \sum_p(\text{Var}(I_p) + \text{Var}(\theta_p)I_p^2) \\\ \sigma_{\Re} = \sqrt{\text{Var}(I_{\Re})} &&\approx && \sqrt{\frac{1}{9} \sum_p(\text{Var}(I_p) \cos^2(\theta_p + \frac{2 \pi}{3}p) + \text{Var}(\theta_p)I_p^2 \sin^2(\theta_p + \frac{2\pi}{3}p))} \\\ \sigma_{\Im} = \sqrt{\text{Var}(I_{\Im})} &&\approx && \sqrt{\frac{1}{9} \sum_p(\text{Var}(I_p) \sin^2(\theta_p + \frac{2 \pi}{3}p) + \text{Var}(\theta_p)I_p^2 \cos^2(\theta_p + \frac{2\pi}{3}p))} \\\ \sigma_{\underline{I}} = \sqrt{\sigma_{\Re}^2 + \sigma_{\Im}^2} &&\approx&& \sqrt{\frac{1}{9} \sum_p(\text{Var}(I_p) + \text{Var}(\theta_p)I_p^2)} \\\ \end{eqnarray}$$

Under the assumption mentioned above that all phases have equal angle and magnitude variance, this reduced to the following.

$$\begin{eqnarray} \sigma_{\Re} &&\approx && \sqrt{\frac{\frac{1}{3} \text{Var}(I_m)\sum_p \cos^2(\theta_p + \frac{2 \pi}{3}p) + \frac{1}{3}\text{Var}(\theta_m)\sum_p(I_p^2 \sin^2(\theta_p + \frac{2\pi}{3}p))}{3}} \\\ \sigma_{\Im} &&\approx && \sqrt{\frac{\frac{1}{3} \text{Var}(I_m) \sum_p \sin^2(\theta_p + \frac{2 \pi}{3}p) + \frac{1}{3}\text{Var}(\theta_m)\sum_p (I_p^2 \cos^2(\theta_p + \frac{2\pi}{3}p))}{3}} \\\ \sigma_{\underline{I}} = \sqrt{\sigma_{\Re}^2 + \sigma_{\Im}^2} &&\approx && \sqrt{\frac{\text{Var}(I_m) + \frac{1}{3}\text{Var}(\theta_m) \sum_p I_p^2}{3}} \end{eqnarray}$$

The outer $$\frac{1}{3}$$ gives rise to the typical $$\frac{1}{\sqrt{3}}$$ often encountered in statistics when aggregating values.

Local angle measurements

For local angle measurements, the angle was obtained asymmetrically and, therefore, the phase shift to obtain the positive sequence has already been accounted for. As a result, the measured values of the different phases can be interpreted as multiple measurements of the same value. Therefore, $$\underline{I}$$ can be interpreted as the mean value of the different measurements:

$$\begin{eqnarray} \underline{I} \coloneqq \underline{I}_1 &&=&& \frac{\sum_p I_p e^{i\theta_p}}{3} \\\ \text{Var}(I_{\Re}) &&\approx && \frac{1}{9} \sum_p(\text{Var}(I_p) \cos^2 \theta_p + \text{Var}(\theta_p)I_p^2 \sin^2\theta_p) \\\ \text{Var}(I_{\Im}) &&\approx && \frac{1}{9} \sum_p(\text{Var}(I_p) \sin^2 \theta_p + \text{Var}(\theta_p)I_p^2 \cos^2 \theta_p) \\\ \text{Var}(\underline{I}) &&\approx && \frac{1}{9} \sum_p(\text{Var}(I_p) + \text{Var}(\theta_p)I_p^2) \\\ \sigma_{\Re} = \sqrt{\text{Var}(I_{\Re})} &&\approx && \sqrt{\frac{1}{9} \sum_p(\text{Var}(I_p) \cos^2 \theta_p + \text{Var}(\theta_p)I_p^2 \sin^2 \theta_p)} \\\ \sigma_{\Im} = \sqrt{\text{Var}(I_{\Im})} &&\approx && \sqrt{\frac{1}{9} \sum_p(\text{Var}(I_p) \sin^2\theta_p + \text{Var}(\theta_p)I_p^2 \cos^2\theta_p)} \\\ \sigma_{\underline{I}} = \sqrt{\sigma_{\Re}^2 + \sigma_{\Im}^2} &&\approx&& \sqrt{\frac{1}{9} \sum_p(\text{Var}(I_p) + \text{Var}(\theta_p)I_p^2)} \\\ \end{eqnarray}$$

Under the assumption mentioned above that all phases have equal angle and magnitude variance, this reduced to the following.

$$\begin{eqnarray} \sigma_{\Re} &&\approx && \sqrt{\frac{\frac{1}{3} \text{Var}(I_m)\sum_p \cos^2\theta_p + \frac{1}{3}\text{Var}(\theta_m)\sum_p(I_p^2 \sin^2\theta_p)}{3}} \\\ \sigma_{\Im} &&\approx && \sqrt{\frac{\frac{1}{3} \text{Var}(I_m) \sum_p \sin^2\theta_p + \frac{1}{3}\text{Var}(\theta_m)\sum_p (I_p^2 \cos^2\theta_p)}{3}} \\\ \sigma_{\underline{I}} = \sqrt{\sigma_{\Re}^2 + \sigma_{\Im}^2} &&\approx && \sqrt{\frac{\text{Var}(I_m) + \frac{1}{3}\text{Var}(\theta_m) \sum_p I_p^2}{3}} \end{eqnarray}$$

The value and variances for symmetric calculations for asymmetric current sensors are therefore the component-wise (phase-wise) aggregation of those for asymmetric calculations (see below).

Asymmetric calculations

$$\begin{eqnarray} && && \underline{I}_p &&\coloneqq && I_p e^{i \theta_p} \\\ I_{\Re,p} &&\coloneqq && \Re\{\underline{I}_p\} &&=&& I_p \cos \theta_p \\\ I_{\Im,p} &&\coloneqq && \Im\{\underline{I}_p\} &&=&& I_p \sin \theta_p \\\ && && \text{Var}(I_{\Re,p}) &&\approx && \text{Var}(I_p) \cos^2 \theta_p + \text{Var}(\theta_p)I_p^2 \sin^2\theta_p \\\ && && \text{Var}(I_{\Im,p}) &&\approx && \text{Var}(I_p) \sin^2 \theta_p + \text{Var}(\theta_p)I_p^2 \cos^2 \theta_p \\\ && && \text{Var}(\underline{I}_p) &&\approx && \text{Var}(I_p) + \text{Var}(\theta_p)I_p^2 \\\ \sigma_{\Re,p} &&=&& \sqrt{\text{Var}(I_{\Re,p})} && \approx && \sqrt{\text{Var}(I_p) \cos^2 \theta_p + \text{Var}(\theta_p)I_p^2 \sin^2\theta_p} \\\ \sigma_{\Im,p} &&=&& \sqrt{\text{Var}(I_{\Im,p})} && \approx && \sqrt{\text{Var}(I_p) \sin^2 \theta_p + \text{Var}(\theta_p)I_p^2 \cos^2 \theta_p} \\\ \sigma_{p} &&=&& \sqrt{\text{Var}(I_p e^{i\theta_p})} \\\ &&=&& \sqrt{\sigma_{\Im,p}^2 + \sigma_{\Re,p}^2} && \approx && \sqrt{\text{Var}(I_p) + \text{Var}(\theta_p)I_p^2} \end{eqnarray}$$

Under the assumption mentioned above that all phases have equal angle and magnitude variance, this reduced to the following.

$$\begin{eqnarray} \sigma_{\Re,p} && \approx && \sqrt{\text{Var}(I_m) \cos^2 \theta_p + \text{Var}(\theta_m)I_p^2 \sin^2\theta_p} \\\ \sigma_{\Im,p} && \approx && \sqrt{\text{Var}(I_m) \sin^2 \theta_p + \text{Var}(\theta_m)I_p^2 \cos^2 \theta_p} \\\ \sigma_{p} = \sqrt{\sigma_{\Im,p}^2 + \sigma_{\Re,p}^2} && \approx && \sqrt{\text{Var}(I_m) + \text{Var}(\theta_m)I_p^2} \\\ \end{eqnarray}$$

[petersalemink95]

@TonyXiang8787 you mention that a sym/asym current_sensor should have i_p_measured and i_q_measured, which means the decoupling lies at the end-user. Wouldn't we in this case want to have i_measured and i_angle_measured, i.e. the current magnitude and angle, where we would do the decoupling? From a user perspective it would be beneficial if they can just provide the measured values, instead of calculating i_p_measured and i_q_measured themselves

[TonyXiang8787]

@TonyXiang8787 you mention that a sym/asym current_sensor should have i_p_measured and i_q_measured, which means the decoupling lies at the end-user. Wouldn't we in this case want to have i_measured and i_angle_measured, i.e. the current magnitude and angle, where we would do the decoupling? From a user perspective it would be beneficial if they can just provide the measured values, instead of calculating i_p_measured and i_q_measured themselves

I agree with you and will adjust the issue. Also we might want to support current phasors, and this interface can be re-used.

[mgovers]

Cfr. offline discussion: the above workout for the current sensor statistics is approved