diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp
index 247a9b676..55fcd1ff3 100644
--- a/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp
+++ b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp
@@ -208,6 +208,51 @@ template <symmetry_tag sym_type> struct PolarComplexRandVar {
 };
 
 namespace statistics {
+template <symmetry_tag sym, template <symmetry_tag> typename RandVarType>
+    requires is_in_list_c<RandVarType<sym>, UniformRealRandVar<sym>, IndependentRealRandVar<sym>,
+                          UniformComplexRandVar<sym>, IndependentComplexRandVar<sym>>
+inline auto scale(RandVarType<sym> const& var, double scale_factor) {
+    return RandVarType<sym>{.value = var.value * scale_factor, .variance = var.variance * abs2(scale_factor)};
+}
+
+template <typename RandVarType>
+    requires is_in_list_c<RandVarType, IndependentRealRandVar<asymmetric_t>, IndependentComplexRandVar<asymmetric_t>>
+inline auto scale(RandVarType const& var, RealValue<asymmetric_t> const& scale_factor) {
+    return RandVarType{.value = var.value * scale_factor, .variance = var.variance * abs2(scale_factor)};
+}
+
+template <symmetry_tag sym, template <symmetry_tag> typename RandVarType>
+    requires is_in_list_c<RandVarType<sym>, UniformComplexRandVar<sym>, IndependentComplexRandVar<sym>>
+inline auto scale(RandVarType<sym> const& var, DoubleComplex const& scale_factor) {
+    return RandVarType<sym>{.value = var.value * scale_factor, .variance = var.variance * abs2(scale_factor)};
+}
+
+inline auto scale(IndependentComplexRandVar<asymmetric_t> const& var, ComplexValue<asymmetric_t> const& scale_factor) {
+    return IndependentComplexRandVar<asymmetric_t>{.value = var.value * scale_factor,
+                                                   .variance = var.variance * abs2(scale_factor)};
+}
+
+template <symmetry_tag sym, typename ScaleType>
+    requires is_in_list_c<ScaleType, RealValue<symmetric_t>, RealValue<asymmetric_t>>
+inline auto scale(DecomposedComplexRandVar<sym> const& var, ScaleType const& scale_factor) {
+    return DecomposedComplexRandVar<sym>{.real_component = scale(var.real_component, scale_factor),
+                                         .imag_component = scale(var.imag_component, scale_factor)};
+}
+
+template <symmetry_tag sym, typename ScaleType>
+    requires(std::same_as<ScaleType, ComplexValue<symmetric_t>> ||
+             (is_asymmetric_v<sym> && std::same_as<ScaleType, ComplexValue<asymmetric_t>>))
+inline auto scale(DecomposedComplexRandVar<sym> const& var, ScaleType const& scale_factor) {
+    ComplexValue<sym> const scaled_value = var.value() * scale_factor;
+    return DecomposedComplexRandVar<sym>{
+        .real_component = {.value = real(scaled_value),
+                           .variance = var.real_component.variance * abs2(real(scale_factor)) +
+                                       var.imag_component.variance * abs2(imag(scale_factor))},
+        .imag_component = {.value = imag(scaled_value),
+                           .variance = var.real_component.variance * abs2(imag(scale_factor)) +
+                                       var.imag_component.variance * abs2(real(scale_factor))}};
+}
+
 // combine multiple random variables of one quantity using Kalman filter
 template <std::ranges::view RandVarsView>
     requires std::same_as<std::ranges::range_value_t<RandVarsView>,
diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp
index c10d31c1e..213fb4aae 100644
--- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp
+++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp
@@ -98,12 +98,6 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
         auto const observability_result =
             necessary_observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure());
 
-        // prepare matrix
-        sub_timer = Timer(calculation_info, 2222, "Prepare matrix, including pre-factorization");
-        prepare_matrix(y_bus, measured_values);
-        // prefactorize
-        sparse_solver_.prefactorize(data_gain_, perm_, observability_result.use_perturbation());
-
         // initialize voltage with initial angle
         sub_timer = Timer(calculation_info, 2223, "Initialize voltages");
         RealValue<sym> const mean_angle_shift = measured_values.mean_angle_shift();
@@ -117,8 +111,21 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
             if (num_iter++ == max_iter) {
                 throw IterationDiverge{max_iter, max_dev, err_tol};
             }
+
+            // prepare matrix
+            sub_timer = Timer(calculation_info, 2222, "Linearize measurements");
+            // get generated (measured/estimated) voltage phasor
+            // with current result voltage angle
+            ComplexValueVector<sym> u = linearize_measurements(output.u, measured_values);
+
+            // prepare matrix
+            sub_timer = Timer(calculation_info, 2222, "Prepare matrix, including pre-factorization");
+            prepare_matrix(y_bus, measured_values, u);
+            // prefactorize
+            sparse_solver_.prefactorize(data_gain_, perm_, observability_result.use_perturbation());
+
             sub_timer = Timer(calculation_info, 2224, "Calculate rhs");
-            prepare_rhs(y_bus, measured_values, output.u);
+            prepare_rhs(y_bus, measured_values, u);
             // solve with prefactorization
             sub_timer = Timer(calculation_info, 2225, "Solve sparse linear equation (pre-factorized)");
             sparse_solver_.solve_with_prefactorized_matrix(data_gain_, perm_, x_rhs_, x_rhs_);
@@ -163,10 +170,12 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
         return ComplexDiagonalTensor<sym>{static_cast<ComplexValue<sym>>(RealValue<sym>{1.0} / value)};
     }
 
-    void prepare_matrix(YBus<sym> const& y_bus, MeasuredValues<sym> const& measured_value) {
+    void prepare_matrix(YBus<sym> const& y_bus, MeasuredValues<sym> const& measured_value,
+                        ComplexValueVector<sym> const& u) {
         MathModelParam<sym> const& param = y_bus.math_model_param();
         IdxVector const& row_indptr = y_bus.row_indptr_lu();
         IdxVector const& col_indices = y_bus.col_indices_lu();
+        std::vector<BranchIdx> const& branch_bus_idx = y_bus.math_topology().branch_bus_idx;
 
         // loop data index, all rows and columns
         for (Idx row = 0; row != n_bus_; ++row) {
@@ -197,9 +206,12 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
                         if (measured_value.has_shunt(obj)) {
                             // G += (-Ys)^H * (variance^-1) * (-Ys)
                             auto const& shunt_power = measured_value.shunt_power(obj);
+                            auto const shunt_current_conj =
+                                statistics::scale(shunt_power, ComplexValue<sym>{1.0 / u[row]});
                             block.g() +=
                                 dot(hermitian_transpose(param.shunt_param[obj]),
-                                    diagonal_inverse(static_cast<IndependentComplexRandVar<sym>>(shunt_power).variance),
+                                    diagonal_inverse(
+                                        static_cast<IndependentComplexRandVar<sym>>(shunt_current_conj).variance),
                                     param.shunt_param[obj]);
                         }
                     }
@@ -214,9 +226,15 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
                             if (std::invoke(has_branch_power_[measured_side], measured_value, obj)) {
                                 // G += Y{side, b0}^H * (variance^-1) * Y{side, b1}
                                 auto const& power = std::invoke(branch_power_[measured_side], measured_value, obj);
+                                // the current needs to be calculated with the voltage of the measured bus side
+                                // NOTE: not the current bus!
+                                Idx const measured_bus = branch_bus_idx[obj][measured_side];
+                                auto const current_conj =
+                                    statistics::scale(power, ComplexValue<sym>{1.0 / u[measured_bus]});
                                 block.g() +=
                                     dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b0]),
-                                        diagonal_inverse(static_cast<IndependentComplexRandVar<sym>>(power).variance),
+                                        diagonal_inverse(
+                                            static_cast<IndependentComplexRandVar<sym>>(current_conj).variance),
                                         param.branch_param[obj].value[measured_side * 2 + b1]);
                             }
                         }
@@ -232,7 +250,9 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
                         // assign variance to diagonal of 3x3 tensor, for asym
                         auto const& injection = measured_value.bus_injection(row);
                         block.r() = ComplexTensor<sym>{static_cast<ComplexValue<sym>>(
-                            -(static_cast<IndependentComplexRandVar<sym>>(injection).variance))};
+                            -(static_cast<IndependentComplexRandVar<sym>>(
+                                  statistics::scale(injection, ComplexValue<sym>{1.0 / u[row]}))
+                                  .variance))};
                     }
                 }
                 // injection measurement not exist
@@ -260,12 +280,9 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
     }
 
     void prepare_rhs(YBus<sym> const& y_bus, MeasuredValues<sym> const& measured_value,
-                     ComplexValueVector<sym> const& current_u) {
+                     ComplexValueVector<sym> const& u) {
         MathModelParam<sym> const& param = y_bus.math_model_param();
-        std::vector<BranchIdx> const branch_bus_idx = y_bus.math_topology().branch_bus_idx;
-        // get generated (measured/estimated) voltage phasor
-        // with current result voltage angle
-        ComplexValueVector<sym> u = linearize_measurements(current_u, measured_value);
+        std::vector<BranchIdx> const& branch_bus_idx = y_bus.math_topology().branch_bus_idx;
 
         // loop all bus to fill rhs
         for (Idx bus = 0; bus != n_bus_; ++bus) {
@@ -286,12 +303,14 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
                 // shunt
                 if (type == YBusElementType::shunt) {
                     if (measured_value.has_shunt(obj)) {
-                        PowerSensorCalcParam<sym> const& m = measured_value.shunt_power(obj);
+                        PowerSensorCalcParam<sym> const& power = measured_value.shunt_power(obj);
+                        auto const current_conj = statistics::scale(power, ComplexValue<sym>{1.0 / u[bus]});
+                        auto const rot_invariant_current_conj =
+                            static_cast<IndependentComplexRandVar<sym>>(current_conj);
                         // eta += (-Ys)^H * (variance^-1) * i_shunt
-                        rhs_block.eta() -=
-                            dot(hermitian_transpose(param.shunt_param[obj]),
-                                diagonal_inverse(static_cast<IndependentComplexRandVar<sym>>(m).variance),
-                                conj(m.value() / u[bus]));
+                        rhs_block.eta() -= dot(hermitian_transpose(param.shunt_param[obj]),
+                                               diagonal_inverse(rot_invariant_current_conj.variance),
+                                               conj(rot_invariant_current_conj.value));
                     }
                 }
                 // branch
@@ -303,23 +322,28 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
                     for (IntS const measured_side : std::array<IntS, 2>{0, 1}) {
                         // has measurement
                         if (std::invoke(has_branch_power_[measured_side], measured_value, obj)) {
-                            PowerSensorCalcParam<sym> const& m =
+                            PowerSensorCalcParam<sym> const& power =
                                 std::invoke(branch_power_[measured_side], measured_value, obj);
                             // the current needs to be calculated with the voltage of the measured bus side
                             // NOTE: not the current bus!
                             Idx const measured_bus = branch_bus_idx[obj][measured_side];
+                            auto const current_conj =
+                                statistics::scale(power, ComplexValue<sym>{1.0 / u[measured_bus]});
+                            auto const rot_invariant_current_conj =
+                                static_cast<IndependentComplexRandVar<sym>>(current_conj);
                             // eta += Y{side, b}^H * (variance^-1) * i_branch_{f, t}
                             rhs_block.eta() +=
                                 dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b]),
-                                    diagonal_inverse(static_cast<IndependentComplexRandVar<sym>>(m).variance),
-                                    conj(m.value() / u[measured_bus]));
+                                    diagonal_inverse(rot_invariant_current_conj.variance),
+                                    conj(rot_invariant_current_conj.value));
                         }
                     }
                 }
             }
             // fill block with injection measurement, need to convert to current
             if (measured_value.has_bus_injection(bus)) {
-                rhs_block.tau() = conj(measured_value.bus_injection(bus).value() / u[bus]);
+                rhs_block.tau() =
+                    conj(statistics::scale(measured_value.bus_injection(bus), ComplexValue<sym>{1.0 / u[bus]}).value());
             }
         }
     }
diff --git a/tests/cpp_unit_tests/test_statistics.cpp b/tests/cpp_unit_tests/test_statistics.cpp
index 4df368497..504433959 100644
--- a/tests/cpp_unit_tests/test_statistics.cpp
+++ b/tests/cpp_unit_tests/test_statistics.cpp
@@ -815,6 +815,249 @@ TEST_CASE("Test statistics") {
     }
 }
 
+TEST_CASE("Test statistics - scale") {
+    using statistics::scale;
+
+    SUBCASE("UniformRealRandVar<symmetric_t> | IndependentRealRandVar<symmetric_t>") {
+        auto const check = []<typename RandVar> {
+            RandVar const var{.value = 1.0, .variance = 2.0};
+            auto const scaled = scale(var, 3.0);
+            CHECK(scaled.value == doctest::Approx(3.0));
+            CHECK(scaled.variance == doctest::Approx(18.0));
+        };
+
+        SUBCASE("UniformRealRandVar<symmetric_t>") { check.template operator()<UniformRealRandVar<symmetric_t>>(); }
+        SUBCASE("IndependentRealRandVar<symmetric_t>") {
+            check.template operator()<IndependentRealRandVar<symmetric_t>>();
+        }
+    }
+    SUBCASE("UniformRealRandVar<asymmetric_t>") {
+        UniformRealRandVar<asymmetric_t> const var{.value = {1.0, 2.0, 3.0}, .variance = 2.0};
+        SUBCASE("Scalar scale") {
+            auto const scaled = scale(var, 3.0);
+            CHECK(scaled.value(0) == doctest::Approx(3.0));
+            CHECK(scaled.value(1) == doctest::Approx(6.0));
+            CHECK(scaled.value(2) == doctest::Approx(9.0));
+            CHECK(scaled.variance == doctest::Approx(18.0));
+        }
+        // scaling asymmetrically would promote the UniformRealRandVar to an IndependentRealRandVar, because the
+        // individual phases scale differently
+    }
+    SUBCASE("IndependentRealRandVar<asymmetric_t>") {
+        auto const var = IndependentRealRandVar<asymmetric_t>{.value = {1.0, 2.0, 3.0}, .variance = {2.0, 3.0, 4.0}};
+        SUBCASE("Scalar scale") {
+            auto const scaled = scale(var, 3.0);
+            CHECK(scaled.value(0) == doctest::Approx(3.0));
+            CHECK(scaled.value(1) == doctest::Approx(6.0));
+            CHECK(scaled.value(2) == doctest::Approx(9.0));
+            CHECK(scaled.variance(0) == doctest::Approx(18.0));
+            CHECK(scaled.variance(1) == doctest::Approx(27.0));
+            CHECK(scaled.variance(2) == doctest::Approx(36.0));
+        }
+        SUBCASE("Asymmetric scale") {
+            auto const scaled = scale(var, RealValue<asymmetric_t>{1.0, 2.0, 3.0});
+            CHECK(scaled.value(0) == doctest::Approx(1.0));
+            CHECK(scaled.value(1) == doctest::Approx(4.0));
+            CHECK(scaled.value(2) == doctest::Approx(9.0));
+            CHECK(scaled.variance(0) == doctest::Approx(2.0));
+            CHECK(scaled.variance(1) == doctest::Approx(12.0));
+            CHECK(scaled.variance(2) == doctest::Approx(36.0));
+        }
+    }
+    SUBCASE("UniformComplexRandVar<symmetric_t> | IndependentComplexRandVar<symmetric_t>") {
+        auto const check = []<typename RandVar> {
+            RandVar const var{.value = {1.0, 2.0}, .variance = 2.0};
+
+            SUBCASE("Real scalar scale") {
+                auto const scaled = scale(var, 3.0);
+                CHECK(real(scaled.value) == doctest::Approx(3.0));
+                CHECK(imag(scaled.value) == doctest::Approx(6.0));
+                CHECK(scaled.variance == doctest::Approx(18.0));
+            }
+            SUBCASE("Complex scalar scale") {
+                auto const scaled = scale(var, DoubleComplex{3.0, 4.0});
+                CHECK(real(scaled.value) == doctest::Approx(3.0 * 1.0 - 4.0 * 2.0));
+                CHECK(imag(scaled.value) == doctest::Approx(3.0 * 2.0 + 4.0 * 1.0));
+                CHECK(scaled.variance == doctest::Approx(2.0 * (3.0 * 3.0 + 4.0 * 4.0)));
+            }
+        };
+        SUBCASE("UniformComplexRandVar<symmetric_t>") {
+            check.template operator()<UniformComplexRandVar<symmetric_t>>();
+        }
+        SUBCASE("IndependentComplexRandVar<symmetric_t>") {
+            check.template operator()<IndependentComplexRandVar<symmetric_t>>();
+        }
+    }
+    SUBCASE("UniformComplexRandVar<asymmetric_t>") {
+        UniformComplexRandVar<asymmetric_t> const var{
+            .value = {RealValue<asymmetric_t>{1.0, 2.0, 3.0}, RealValue<asymmetric_t>{4.0, 5.0, 6.0}}, .variance = 2.0};
+        SUBCASE("Real scalar scale") {
+            auto const scaled = scale(var, 3.0);
+            CHECK(real(scaled.value(0)) == doctest::Approx(3.0));
+            CHECK(real(scaled.value(1)) == doctest::Approx(6.0));
+            CHECK(real(scaled.value(2)) == doctest::Approx(9.0));
+            CHECK(imag(scaled.value(0)) == doctest::Approx(12.0));
+            CHECK(imag(scaled.value(1)) == doctest::Approx(15.0));
+            CHECK(imag(scaled.value(2)) == doctest::Approx(18.0));
+            CHECK(scaled.variance == doctest::Approx(18.0));
+        }
+        SUBCASE("Complex scalar scale") {
+            DoubleComplex const scale_factor{3.0, 4.0};
+            auto const scaled = scale(var, scale_factor);
+            for (auto i = 0; i < 3; ++i) {
+                CHECK(real(scaled.value(i)) ==
+                      real(scale(UniformComplexRandVar<symmetric_t>{.value = var.value(i), .variance = var.variance},
+                                 scale_factor)
+                               .value));
+                CHECK(real(scaled.value(i)) ==
+                      real(scale(UniformComplexRandVar<symmetric_t>{.value = var.value(i), .variance = var.variance},
+                                 scale_factor)
+                               .value));
+            }
+            CHECK(scaled.variance == doctest::Approx(2.0 * (3.0 * 3.0 + 4.0 * 4.0)));
+        }
+        // scaling asymmetrically would promote the UniformComplexRandVar to an IndependentComplexRandVar, because the
+        // individual phases scale differently
+    }
+    SUBCASE("IndependentComplexRandVar<asymmetric_t>") {
+        auto const var = IndependentComplexRandVar<asymmetric_t>{
+            .value = {RealValue<asymmetric_t>{1.0, 2.0, 3.0}, RealValue<asymmetric_t>{4.0, 5.0, 6.0}},
+            .variance = {2.0, 3.0, 4.0}};
+        auto const individual_phases = [&var] {
+            std::array<UniformComplexRandVar<symmetric_t>, 3> result;
+            for (auto i = 0; i < 3; ++i) {
+                result[i] = {.value = var.value(i), .variance = var.variance(i)};
+            }
+            return result;
+        }();
+        SUBCASE("Real scalar scale") {
+            auto const scale_factor = 3.0;
+            auto const scaled = scale(var, scale_factor);
+            for (auto i = 0; i < 3; ++i) {
+                auto const expected = scale(individual_phases[i], scale_factor);
+                CHECK(real(scaled.value(i)) == doctest::Approx(real(expected.value)));
+                CHECK(imag(scaled.value(i)) == doctest::Approx(imag(expected.value)));
+                CHECK(scaled.variance(i) == doctest::Approx(expected.variance));
+            }
+        }
+        SUBCASE("Complex scalar scale") {
+            DoubleComplex const scale_factor{3.0, 4.0};
+            auto const scaled = scale(var, scale_factor);
+            for (auto i = 0; i < 3; ++i) {
+                auto const expected = scale(individual_phases[i], scale_factor);
+                CHECK(real(scaled.value(i)) == doctest::Approx(real(expected.value)));
+                CHECK(imag(scaled.value(i)) == doctest::Approx(imag(expected.value)));
+                CHECK(scaled.variance(i) == doctest::Approx(expected.variance));
+            }
+        }
+        SUBCASE("Real asymmetric scale") {
+            RealValue<asymmetric_t> const scale_factor{3.0, 4.0, 5.0};
+            auto const scaled = scale(var, scale_factor);
+            for (auto i = 0; i < 3; ++i) {
+                auto const expected = scale(individual_phases[i], scale_factor[i]);
+                CHECK(real(scaled.value(i)) == doctest::Approx(real(expected.value)));
+                CHECK(imag(scaled.value(i)) == doctest::Approx(imag(expected.value)));
+                CHECK(scaled.variance(i) == doctest::Approx(expected.variance));
+            }
+        }
+        SUBCASE("Complex asymmetric scale") {
+            ComplexValue<asymmetric_t> const scale_factor{RealValue<asymmetric_t>{3.0, 4.0, 5.0},
+                                                          RealValue<asymmetric_t>{6.0, 7.0, 8.0}};
+            auto const scaled = scale(var, scale_factor);
+            for (auto i = 0; i < 3; ++i) {
+                auto const expected = scale(individual_phases[i], scale_factor(i));
+                CHECK(real(scaled.value(i)) == doctest::Approx(real(expected.value)));
+                CHECK(imag(scaled.value(i)) == doctest::Approx(imag(expected.value)));
+                CHECK(scaled.variance(i) == doctest::Approx(expected.variance));
+            }
+        }
+    }
+    SUBCASE("DecomposedComplexRandVar<symmetric_t>") {
+        DecomposedComplexRandVar<symmetric_t> const var{.real_component = {.value = 1.0, .variance = 2.0},
+                                                        .imag_component = {.value = 4.0, .variance = 5.0}};
+
+        SUBCASE("Real scalar scale") {
+            auto const scale_factor = 3.0;
+            auto const scaled = scale(var, scale_factor);
+            CHECK(scaled.real_component.value == scale(var.real_component, scale_factor).value);
+            CHECK(scaled.imag_component.value == scale(var.imag_component, scale_factor).value);
+            CHECK(scaled.real_component.variance == scale(var.real_component, scale_factor).variance);
+            CHECK(scaled.imag_component.variance == scale(var.imag_component, scale_factor).variance);
+        }
+        SUBCASE("Complex scalar scale") {
+            auto const scale_factor = DoubleComplex{3.0, 4.0};
+            auto const scaled = scale(var, scale_factor);
+            CHECK(scaled.real_component.value == real(var.value() * scale_factor));
+            CHECK(scaled.imag_component.value == imag(var.value() * scale_factor));
+            CHECK(scaled.real_component.variance ==
+                  real(scale_factor) * real(scale_factor) * var.real_component.variance +
+                      imag(scale_factor) * imag(scale_factor) * var.imag_component.variance);
+            CHECK(scaled.imag_component.variance ==
+                  real(scale_factor) * real(scale_factor) * var.imag_component.variance +
+                      imag(scale_factor) * imag(scale_factor) * var.real_component.variance);
+        }
+    }
+    SUBCASE("DecomposedComplexRandVar<asymmetric_t>") {
+        DecomposedComplexRandVar<asymmetric_t> const var{
+            .real_component = {.value = {1.0, 2.0, 3.0}, .variance = {2.0, 3.0, 4.0}},
+            .imag_component = {.value = {4.0, 5.0, 6.0}, .variance = {5.0, 6.0, 7.0}}};
+
+        SUBCASE("Real scalar scale") {
+            auto const scale_factor = 3.0;
+            auto const scaled = scale(var, scale_factor);
+
+            for (int i = 0; i < 3; ++i) {
+                CHECK(scaled.real_component.value(i) == scale(var.real_component, scale_factor).value(i));
+                CHECK(scaled.imag_component.value(i) == scale(var.imag_component, scale_factor).value(i));
+                CHECK(scaled.real_component.variance(i) == scale(var.real_component, scale_factor).variance(i));
+                CHECK(scaled.imag_component.variance(i) == scale(var.imag_component, scale_factor).variance(i));
+            }
+        }
+        SUBCASE("Real asymmetric scale") {
+            auto const scale_factor = RealValue<asymmetric_t>{3.0, 4.0, 5.0};
+            auto const scaled = scale(var, scale_factor);
+
+            for (int i = 0; i < 3; ++i) {
+                CHECK(scaled.real_component.value(i) == scale(var.real_component, scale_factor).value(i));
+                CHECK(scaled.imag_component.value(i) == scale(var.imag_component, scale_factor).value(i));
+                CHECK(scaled.real_component.variance(i) == scale(var.real_component, scale_factor).variance(i));
+                CHECK(scaled.imag_component.variance(i) == scale(var.imag_component, scale_factor).variance(i));
+            }
+        }
+        SUBCASE("Complex scalar scale") {
+            auto const scale_factor = DoubleComplex{3.0, 4.0};
+            auto const scaled = scale(var, scale_factor);
+
+            for (int i = 0; i < 3; ++i) {
+                CHECK(scaled.real_component.value(i) == real(var.value()(i) * scale_factor));
+                CHECK(scaled.imag_component.value(i) == imag(var.value()(i) * scale_factor));
+                CHECK(scaled.real_component.variance(i) ==
+                      real(scale_factor) * real(scale_factor) * var.real_component.variance(i) +
+                          imag(scale_factor) * imag(scale_factor) * var.imag_component.variance(i));
+                CHECK(scaled.imag_component.variance(i) ==
+                      real(scale_factor) * real(scale_factor) * var.imag_component.variance(i) +
+                          imag(scale_factor) * imag(scale_factor) * var.real_component.variance(i));
+            }
+        }
+        SUBCASE("Complex asymmetric scale") {
+            auto const scale_factor = ComplexValue<asymmetric_t>{RealValue<asymmetric_t>{3.0, 4.0, 5.0},
+                                                                 RealValue<asymmetric_t>{6.0, 7.0, 8.0}};
+            auto const scaled = scale(var, scale_factor);
+
+            for (int i = 0; i < 3; ++i) {
+                CHECK(scaled.real_component.value(i) == real(var.value()(i) * scale_factor(i)));
+                CHECK(scaled.imag_component.value(i) == imag(var.value()(i) * scale_factor(i)));
+                CHECK(scaled.real_component.variance(i) ==
+                      real(scale_factor(i)) * real(scale_factor(i)) * var.real_component.variance(i) +
+                          imag(scale_factor(i)) * imag(scale_factor(i)) * var.imag_component.variance(i));
+                CHECK(scaled.imag_component.variance(i) ==
+                      real(scale_factor(i)) * real(scale_factor(i)) * var.imag_component.variance(i) +
+                          imag(scale_factor(i)) * imag(scale_factor(i)) * var.real_component.variance(i));
+            }
+        }
+    }
+}
+
 TEST_CASE("Test statistics - combine") {
     using statistics::combine;
     using std::views::take;
diff --git a/tests/data/state_estimation/line-power-sensor/params.json b/tests/data/state_estimation/line-power-sensor/params.json
index 47b5ca68b..fb888eb3a 100644
--- a/tests/data/state_estimation/line-power-sensor/params.json
+++ b/tests/data/state_estimation/line-power-sensor/params.json
@@ -3,6 +3,6 @@
   "rtol": 1e-03,
   "atol": {
     "default": 1e-05,
-    ".+_residual": 0.1
+    ".+_residual": 0.2
   }
 }
\ No newline at end of file
diff --git a/tests/data/state_estimation/single-line-load-il/sym_output.json b/tests/data/state_estimation/single-line-load-il/sym_output.json
index 7e34500c8..c0f22327a 100644
--- a/tests/data/state_estimation/single-line-load-il/sym_output.json
+++ b/tests/data/state_estimation/single-line-load-il/sym_output.json
@@ -8,8 +8,8 @@
       {
         "id": 21,
         "energized": 1,
-        "p_from": 1501363.636462,
-        "p_to": -1499319.107870
+        "p_from": 1500680.11502645,
+        "p_to": -1498637.44762342
       }
     ],
     "node": [
@@ -28,14 +28,14 @@
       {
         "id": 41,
         "energized": 1,
-        "p": 1501363.636462
+        "p": 1500680.11502645
       }
     ],
     "sym_load": [
       {
         "id": 31,
         "energized": 1,
-        "p": 1499319.107870,
+        "p": 1498637.44762342,
         "q": 0.0
       }
     ]