Documentation: add current sensor + power sensor variance diff to ILSE docs

docs/user_manual/calculations.md

@@ -57,10 +57,13 @@ Output:
 - Power flow through branches
 - Deviation between measurement values and estimated state
 
-In order to perform a state estimation, the system should be observable. If the system is not observable, the calculation will raise either a `NotObservableError` or
-a `SparseMatrixError`.
-In short, meeting the requirement of observability indicates that the system is either an overdetermined system (when the number of measurements is larger than the number of unknowns.
-For each node, there are two unknowns, `u` and `u_angle`. Due to the relative nature of `u_angle` (relevant only in systems with at least two nodes), in total the following conditions should be met:
+In order to perform a state estimation, the system should be observable. If the system is not observable,
+the calculation will raise either a `NotObservableError` or a `SparseMatrixError`.
+In short, meeting the requirement of observability indicates that the system is either an overdetermined
+system (when the number of independent measurements is larger than the number of unknowns) or an exactly
+determined system (the number of independent measurements equals the number of unknowns).
+For each node, there are two unknowns, `u` and `u_angle`. Due to the relative nature of `u_angle`
+(relevant only in systems with at least two nodes), in total the following conditions should be met:
 
 $$
     \begin{eqnarray}
@@ -82,11 +85,15 @@ The number of measurements can be found by taking the sum of the following:
 - two times the number of nodes with a voltage sensor with magnitude and angle
 - two times the number of nodes without appliances connected
 - two times the number of nodes where all connected appliances are measured by a power sensor
-- two times the number of branches with a power sensor
+- two times the number of branches with a power sensor and/or a current sensor
 
 ```{note}
 Having enough measurements does not necessarily mean that the system is observable. The location of the measurements is
-also of importance. Additionally, there should be at least one voltage measurement.
+also of importance, i.e., the measurements should be topologically independent. Additionally, there should be at least one voltage measurement.
+```
+
+```{note}
+Global angle current measurements require at least one voltage angle measurement to make sense.
 ```
 
 ```{warning}
@@ -99,9 +106,9 @@ In observable systems this helps better outputting correct results. On the other
 Based on the requirements of observability mentioned above, users need to satisfy at least the following conditions for state estimation calculation in `power-grid-model`.
 
 - `n_voltage_sensor >= 1`
-- If no voltage phasor sensors are available, then the following conditions should be satisfied:  `n_unique_power_sensor >= n_bus - 1`. Otherwise: `n_unique_power_sensor + n_voltage_sensor_with_phasor >= n_bus`
+- If no voltage phasor sensors are available, then the following conditions should be satisfied:  `n_unique_power_or_current_sensor >= n_bus - 1`. Otherwise: `n_unique_power_or_current_sensor + n_voltage_sensor_with_phasor >= n_bus`
 
-`n_unique_power_sensor` can be calculated as sum of following:
+`n_unique_power_or_current_sensor` can be calculated as sum of following:
 
 - Zero injection or zero power flow constraint if present for all nodes.
 - Complete injections for all nodes: All appliances in a node are measured or a node injection sensor is present. Either of them counts as one.
@@ -466,8 +473,35 @@ $$
     \end{eqnarray}
 $$
 
+- Branch current with global angle: Linear WLS requires a complex current phasor. The global angle current measurement captures
+the phase offset relative to the same predetermined reference phase against which the voltage angle is measured. As a result,
+at least one voltage measurement with angle measurement is required in the grid if a global angle current sensor is used.
+
+$$
+    \begin{eqnarray}
+            \underline{I} = I_i \cdot e^{j \theta_i}
+    \end{eqnarray}
+$$
+
+- Branch current with local angle: Linear WLS requires a complex current phasor. Sometimes, (accurate) voltage measurements are
+not available for a branch, which means no power measurement is possible. Using only the current amplitude is difficult.
+However, it may be possible to obtain a reasonably accurately measurement of the current amplitude and the relative phase angle
+to the voltage. In this way, we still have enough information for the state estimation. The resulting local current phasor
+$\underline{I}_{\text{local}}$ can be translated to the global current phasor (see above) for iterative linear state estimation.
+Given the measured (linearized) voltage phasor, the current phasor is calculated as follows:
+
+$$
+    \begin{equation}
+        \underline{I} = \underline{I}_{\text{local}}^{*} \frac{\underline{U}}{|\underline{U}|}
+        = \underline{I}_{\text{local}}^{*} \cdot e^{j \theta}
+    \end{equation}
+$$
+
+where $\underline{U}$ is either measured voltage magnitude at the bus or assumed unity magnitude, with the intrinsic phase shift.
+$\theta$ is the phase angle of the voltage.
+
 - Branch/shunt power flow: Linear WLS requires a complex current phasor. To make this translation, the voltage at the terminal should
-also be measured, otherwise the nominal voltage with zero angle is used as an estimation. With the measured (linearized) voltage
+also be measured, otherwise the nominal voltage with zero angle is used as an estimation. Given the measured (linearized) voltage
 phasor, the current phasor is calculated as follows:
 
 $$
@@ -506,12 +540,28 @@ be the slack bus, which is connected to the external network (source). $\underli
   - Normalize the voltage phasor angle by setting the angle of the slack bus to zero:
   - If the maximum deviation between $\underline{U}^{(k)}$ and $\underline{U}^{(k-1)}$ is smaller than the error tolerance $\epsilon$,
   stop the iteration. Otherwise, continue until the maximum number of iterations is reached.
-  
+
 In the iteration process, the phase angle of voltages at each bus is updated using the last iteration;
 the system error of the phase shift converges to zero. Because the matrix is pre-built and
 pre-factorized, the computation cost of each iteration is much smaller than for the [Newton-Raphson](#newton-raphson-state-estimation)
 method, where the Jacobian matrix needs to be constructed and factorized each time.
 
+Because the matrix depends on the variances of the linearized current measurements, pre-factorization
+requires assuming that the magnitudes of the voltages remain constant throughout the iteration process.
+This assumption is, of course, an approximation, and potentially results an inaccurate relative weighing
+of the power measurements. However, the difference between the solution and the actual grid state is
+usually small, because the error introduced by this is typically dominated by the following contributing
+factors (see also [this thread](https://github.com/PowerGridModel/power-grid-model/pull/951#issuecomment-2805154436)).
+
+- voltages in the distribution grid are usually close to 1 p.u..
+- power sensor errors are usually best-effort estimates.
+- the interpretation of power measurements as current measurements is an approximation.
+
+```{warning}
+In short, the pre-factorization may produce results that may deviate from the actual grid state.
+If higher precision is desired, please consider using [Newton-Raphson state estimation](#newton-raphson-state-estimation) instead.
+```
+
 ```{warning}
 The algorithm will assume angles to be zero by default (see the details about voltage sensors). This produces more correct outputs when the system is observable, but will prevent the calculation from raising an exception, even if it is unobservable, therefore giving faulty results.
 ```
@@ -527,7 +577,7 @@ However, the Jacobian matrix needs to be calculated every iteration and cannot b
 
 The Newton-Raphson method considers all measurements to be independent real measurements.
 I.e., $\sigma_P$, $\sigma_Q$ and $\sigma_U$ are all independent values.
-The rationale behind to calculation is similar to that of the [Newton-Raphson for power flow](#newton-raphson-power-flow).
+The rationale behind to calculation is similar to that of the [Newton-Raphson method for power flow](#newton-raphson-power-flow).
 Consequently, the iteration process differs slightly from that of [iterative linear state estimation](#iterative-linear-state-estimation), as shown below.
 
 - Initialization: let $\boldsymbol{U}^{(k)}$ be the column vector of the estimated voltage magnitude and $\boldsymbol{\theta}^{(k)}$ the column vector of the
@@ -550,6 +600,10 @@ We initialize $\boldsymbol{U}^{(0)}$ and $\boldsymbol{\theta}^{(k)}$ as follows:
 As for the [iterative linear](#iterative-linear-state-estimation) approach, during iterations, phase angles of voltage at each bus are updated using ones from the previous iteration.
 The system error of the phase shift converges to zero.
 
+```{note}
+Newton-Raphson state estimation does not support current measurements at this moment. See also [this issue](https://github.com/PowerGridModel/power-grid-model/issues/765).
+```
+
 ```{warning}
 The algorithm will assume angles to be zero by default (see the details about voltage sensors). In observable systems this helps better outputting correct results. On the other hand with unobservable systems, exceptions raised from calculations due to faulty results will be prevented.
 ```

docs/user_manual/components.md

@@ -721,6 +721,98 @@ $$
    \end{eqnarray}
 $$
 
+### Generic Current Sensor
+
+```{warning}
+At the time of writing, this feature is still experimental and is not yet publicly available.
+```
+
+* type name: `generic_current_sensor`
+
+`current_sensor` is an abstract class for symmetric and asymmetric current sensor and is derived from
+{hoverxreftooltip}`user_manual/components:sensor`. It measures the magnitude and angle of the current flow of a terminal.
+The terminal is connecting the from/to end of a `branch` (except `link`) and a `node`.
+
+```{note}
+Due to the high admittance of a `link` it is chosen that a current sensor cannot be coupled to a `link`, even though a link is a `branch`
+```
+
+##### Input
+
+| name                     | data type                                                                     | unit       | description                                                                                                                               |              required               |  update  |                     valid values                     |
+| ------------------------ | ----------------------------------------------------------------------------- | ---------- | ----------------------------------------------------------------------------------------------------------------------------------------- | :---------------------------------: | :------: | :--------------------------------------------------: |
+| `measured_terminal_type` | {py:class}`MeasuredTerminalType <power_grid_model.enum.MeasuredTerminalType>` | -          | indicate the side of the `branch`                                                                                                         |              &#10004;               | &#10060; | the terminal type should match the `measured_object` |
+| `angle_measurement_type` | {py:class}`AngleMeasurementType <power_grid_model.enum.AngleMeasurementType>` | -          | indicate whether the measured angle is a global angle or a local angle                                                                    |              &#10004;               | &#10060; | the terminal type should match the `measured_object` |
+| `i_sigma`                | `double`                                                                      | ampere (A) | standard deviation of the current (`i`) measurement error. Usually this is the absolute measurement error range divided by 3.             | &#10060; see the explanation below. | &#10004; |                        `> 0`                         |
+| `i_angle_sigma`          | `double`                                                                      | rad        | standard deviation of the current (`i`) phase angle measurement error. Usually this is the absolute measurement error range divided by 3. | &#10060; see the explanation below. | &#10004; |                        `> 0`                         |
+
+#### Current Sensor Concrete Types
+
+```{warning}
+At the time of writing, this feature is still experimental and is not yet publicly available.
+```
+
+There are two concrete types of current sensor. They share similar attributes:
+the meaning of `RealValueInput` is different, as shown in the table below.
+
+| type name             | meaning of `RealValueInput` |
+| --------------------- | --------------------------- |
+| `sym_current_sensor`  | `double`                    |
+| `asym_current_sensor` | `double[3]`                 |
+
+##### Input
+
+| name               | data type        | unit       | description                               |              required              |  update  |
+| ------------------ | ---------------- | ---------- | ----------------------------------------- | :--------------------------------: | :------: |
+| `i_measured`       | `RealValueInput` | ampere (A) | measured current (`i`) magnitude          | &#10024; only for state estimation | &#10004; |
+| `i_angle_measured` | `RealValueInput` | rad        | measured phase angle of the current (`i`) | &#10024; only for state estimation | &#10004; |
+
+See the documentation on [state estimation calculation methods](calculations.md#state-estimation-algorithms) for details per method on how the variances are taken into account for both the global and local angle measurement types and for the individual phases.
+
+##### Steady state output
+
+```{note}
+A sensor only has output for state estimation. For other calculation types, sensor output is undefined.
+```
+
+| name               | data type         | unit       | description                                                                                 |
+| ------------------ | ----------------- | ---------- | ------------------------------------------------------------------------------------------- |
+| `i_residual`       | `RealValueOutput` | ampere (A) | residual value between measured current (`i`) and calculated current (`i`)                  |
+| `i_angle_residual` | `RealValueOutput` | rad        | residual value between measured phase angle and calculated phase angle of the current (`i`) |
+
+#### Electric Model
+
+`Generic Current Sensor` is modeled by following equations:
+
+$$
+   \begin{eqnarray}
+        & p_{\text{residual}} = p_{\text{measured}} - p_{\text{state}} \\
+        & q_{\text{residual}} = q_{\text{measured}} - q_{\text{state}}
+   \end{eqnarray}
+$$
+
+Internally, the current phasor measurement is decomposed into a separate real and imaginary component
+per phase, with each their own variances, which are estimated using the common linearization approximation:
+
+$$
+   \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}
+$$
+
+The following illustrates how this works for `sym_current_sensor`s in symmetric calculations.
+See also [the full mathematical workout](https://github.com/PowerGridModel/power-grid-model/issues/547).
+
+$$
+   \begin{eqnarray}
+        & \real\left\{I\right\}\right = I \cos\theta \\
+        & \imag\left\{I\right\}\right = I \sin\theta \\
+        & \text{Var}\left(\real\left\{I\right\}\right) = \sigma_i^2 \cos^2\theta + I^2 \sigma_{\theta}^2\sin^2\theta \\
+        & \text{Var}\left(\imag\left\{I\right\}\right) = \sigma_i^2 \sin^2\theta + I^2 \sigma_{\theta}^2\cos^2\theta \\
+   \end{eqnarray}
+$$
+
 ## Fault
 
 * type name: `fault`