Jd/fix reduction indexing

src/core/network_model.jl

@@ -129,10 +129,21 @@
         model.PTDF_matrix =
             PNM.PTDF(sys; reduce_radial_branches = model.reduce_radial_branches)
     end
+
+    if !model.reduce_radial_branches && !isempty(model.PTDF_matrix.radial_network_reduction)
+        throw(
+            IS.ConflictingInputsError(
+                "The provided PTDF Matrix has reduced radial branches and mismatches the network \\
+                model specification radial_network_reduction = false. Set the keyword argument \\
+                radial_network_reduction = true in your network model"),
+        )
+    end
+
     if model.reduce_radial_branches
         @assert !isempty(model.PTDF_matrix.radial_network_reduction)
         model.radial_network_reduction = model.PTDF_matrix.radial_network_reduction
     end
+
     get_PTDF_matrix(model).subnetworks
     model.subnetworks = deepcopy(get_PTDF_matrix(model).subnetworks)
     if length(model.subnetworks) > 1

src/devices_models/devices/common/add_to_expression.jl

@@ -236,14 +236,16 @@ function add_to_expression!(
     variable = get_variable(container, U(), V)
     expression = get_expression(container, T(), PSY.ACBus)
     radial_network_reduction = get_radial_network_reduction(network_model)
-    for d in devices, t in get_time_steps(container)
+    for d in devices
         name = PSY.get_name(d)
         bus_no = PNM.get_mapped_bus_number(radial_network_reduction, PSY.get_bus(d))
-        _add_to_jump_expression!(
-            expression[bus_no, t],
-            variable[name, t],
-            get_variable_multiplier(U(), V, W()),
-        )
+        for t in get_time_steps(container)
+            _add_to_jump_expression!(
+                expression[bus_no, t],
+                variable[name, t],
+                get_variable_multiplier(U(), V, W()),
+            )
+        end
     end
     return
 end
@@ -1462,6 +1464,30 @@ function add_to_expression!(
     return
 end
 
+function add_to_expression!(
+    container::OptimizationContainer,
+    ::Type{T},
+    ::Type{U},
+    sys::PSY.System,
+    ::NetworkModel{W},
+) where {
+    T <: ActivePowerBalance,
+    U <: Union{SystemBalanceSlackUp, SystemBalanceSlackDown},
+    W <: AreaBalancePowerModel,
+}
+    variable = get_variable(container, U(), PSY.Area)
+    expression = get_expression(container, T(), PSY.Area)
+    @assert_op length(axes(variable, 1)) == length(axes(expression, 1))
+    for t in get_time_steps(container), n in axes(expression, 1)
+        _add_to_jump_expression!(
+            expression[n, t],
+            variable[n, t],
+            get_variable_multiplier(U(), PSY.Area, W),
+        )
+    end
+    return
+end
+
 function add_to_expression!(
     container::OptimizationContainer,
     ::Type{T},

test/test_network_constructors.jl

@@ -712,6 +712,40 @@ end
     end
 end
 
+@testset "2 Areas AreaBalance PowerModel - with slacks" begin
+    c_sys = build_system(PSITestSystems, "c_sys5_uc")
+    # Extend the system with two areas
+    areas = [Area("Area_1", 0, 0, 0), Area("Area_2", 0, 0, 0)]
+    add_components!(c_sys, areas)
+    for (i, comp) in enumerate(get_components(ACBus, c_sys))
+        (i < 3) ? set_area!(comp, areas[1]) : set_area!(comp, areas[2])
+    end
+    # Deactivate generators on Area 1: as there is no area interchange defined,
+    # slacks will be required for feasibility
+    for gen in get_components(x -> (get_area(get_bus(x)) == areas[1]), Generator, c_sys)
+        set_available!(gen, false)
+    end
+
+    template = get_thermal_dispatch_template_network(
+        NetworkModel(AreaBalancePowerModel; use_slacks = true),
+    )
+    ps_model = DecisionModel(template, c_sys; optimizer = HiGHS_optimizer)
+
+    @test build!(ps_model; output_dir = mktempdir(; cleanup = true)) ==
+          PSI.ModelBuildStatus.BUILT
+    @test solve!(ps_model) == PSI.RunStatus.SUCCESSFULLY_FINALIZED
+
+    opt_container = PSI.get_optimization_container(ps_model)
+    copper_plate_constraints =
+        PSI.get_constraint(opt_container, CopperPlateBalanceConstraint(), PSY.Area)
+    @test size(copper_plate_constraints) == (2, 24)
+
+    results = OptimizationProblemResults(ps_model)
+    slacks_up = read_variable(results, "SystemBalanceSlackUp__Area")
+    @test all(slacks_up[!, "Area_1"] .> 0.0)
+    @test all(slacks_up[!, "Area_2"] .≈ 0.0)
+end
+
 @testset "2 Areas AreaBalance PowerModel" begin
     c_sys = PSB.build_system(PSISystems, "two_area_pjm_DA")
     transform_single_time_series!(c_sys, Hour(24), Hour(1))