diff --git a/src/topupopt/problems/esipp/model.py b/src/topupopt/problems/esipp/model.py
index a378cc8edd15b74a2cea1004c118e276c8cebbf8..508ed3e9ddf2c71b4d1f277de99e8889db205e4c 100644
--- a/src/topupopt/problems/esipp/model.py
+++ b/src/topupopt/problems/esipp/model.py
@@ -957,20 +957,6 @@ def create_model(
         initialize=(init_set_GLLJH_arc_inv_sos1 if enable_initialisation else None),
     )
 
-    # set of GLLJ-indexed GLLJH tuples for new arcs modelled using SOS1
-
-    def init_set_GLLJH_arc_inv_sos1_gllj(m, g, l1, l2, j):
-        return ((g, l1, l2, j, h) for h in m.set_H_gllj[(g, l1, l2, j)])
-
-    model.set_GLLJH_arc_inv_sos1_gllj = pyo.Set(
-        model.set_GLLJ_arc_inv_sos1,
-        dimen=5,
-        within=model.set_GLLJH_arc_inv_sos1,
-        initialize=(
-            init_set_GLLJH_arc_inv_sos1_gllj if enable_initialisation else None
-        ),
-    )
-
     # *************************************************************************
 
     # flow sense determination using SOS1
@@ -1309,18 +1295,6 @@ def create_model(
         dimen=2, within=model.set_TH, initialize=init_set_TH_arc_inv_sos1
     )
 
-    # set of t-indexed TH tuples for groups of arcs relying on SOS1
-
-    def init_set_TH_arc_inv_sos1_t(m, t):
-        return ((t, h) for h in m.set_H_t[t])
-
-    model.set_TH_arc_inv_sos1_t = pyo.Set(
-        model.set_T_sos1,
-        dimen=2,
-        within=model.set_TH_arc_inv_sos1,
-        initialize=init_set_TH_arc_inv_sos1_t,
-    )
-
     # minimum cost of a group of arcs
 
     model.param_c_arc_min_th = pyo.Param(model.set_TH, within=pyo.NonNegativeReals)
@@ -2429,14 +2403,21 @@ def create_model(
 
     # *************************************************************************
 
-    # SOS1 constraints for arc selection
-
+    # SOS1 constraints for arc selection    
+    def rule_constr_arc_sos1(m, g, l1, l2, j):
+        var_list = [
+            m.var_delta_arc_inv_glljh[(g, l1, l2, j, h)]
+            for h in m.set_H_gllj[(g, l1, l2, j)]
+            ]
+        weight_list = [
+            m.param_arc_inv_sos1_weights_glljh[(g, l1, l2, j, h)]
+            for h in m.set_H_gllj[(g, l1, l2, j)]
+            ]
+        return (var_list, weight_list)
     model.constr_arc_sos1 = pyo.SOSConstraint(
-        model.set_GLLJ_arc_inv_sos1,  # for (directed and undirected) new arcs
-        var=model.var_delta_arc_inv_glljh,  # set_GLLJH_sgl indexes the variables
-        index=model.set_GLLJH_arc_inv_sos1_gllj,  # key: GLLJ; value: GLLJH
-        weights=model.param_arc_inv_sos1_weights_glljh,  # key: GLLJH; alue: weight
-        sos=1,
+        model.set_GLLJ_arc_inv_sos1,
+        rule=rule_constr_arc_sos1,
+        sos=1
     )
 
     # *************************************************************************
@@ -2566,16 +2547,29 @@ def create_model(
     # *************************************************************************
 
     # SOS1 constraints for flow sense determination (undirected arcs)
-
+    # model.constr_sns_sos1 = pyo.SOSConstraint(
+    #     model.set_GLLJ_und,  # one constraint per undirected arc
+    #     model.set_QK,  # and time interval
+    #     var=model.var_zeta_sns_glljqk,  # set_GLLJ_und_red and set_K
+    #     index=model.set_GLLJQK_und_sns_sos1_red_gllj,  #
+    #     weights=model.param_arc_sns_sos1_weights_glljqk,
+    #     sos=1,
+    # )
+    # TODO: make the weight list dependent in model options rather than literal
+    def rule_constr_sns_sos1(m, g, l1, l2, j, q, k):
+        var_list = [
+            m.var_zeta_sns_glljqk[(g, l1, l2, j, q, k)],
+            m.var_zeta_sns_glljqk[(g, l2, l1, j, q, k)]
+            ]
+        weight_list = [1, 2] if True else [2, 1]
+        return (var_list, weight_list)
     model.constr_sns_sos1 = pyo.SOSConstraint(
-        model.set_GLLJ_und,  # one constraint per undirected arc
-        model.set_QK,  # and time interval
-        var=model.var_zeta_sns_glljqk,  # set_GLLJ_und_red and set_K
-        index=model.set_GLLJQK_und_sns_sos1_red_gllj,  #
-        weights=model.param_arc_sns_sos1_weights_glljqk,
-        sos=1,
+        model.set_GLLJ_und,
+        model.set_QK,
+        rule=rule_constr_sns_sos1,
+        sos=1
     )
-
+    
     # *************************************************************************
 
     # static losses
@@ -2935,13 +2929,20 @@ def create_model(
     # *************************************************************************
 
     # SOS1 constraints for arc group selection
-
+    def rule_constr_arc_group_sos1(m, t):
+        var_list = [
+            m.var_delta_arc_inv_th[(t, h)]
+            for h in m.set_H_t[t]
+            ]
+        weight_list = [
+            m.param_arc_inv_sos1_weights_th[(t, h)]
+            for h in m.set_H_t[t]
+            ]
+        return (var_list, weight_list)
     model.constr_arc_group_sos1 = pyo.SOSConstraint(
-        model.set_T_sos1,  # for all groups using sos1
-        var=model.var_delta_arc_inv_th,  # set_TH indexes the variables
-        index=model.set_TH_arc_inv_sos1_t,  # key: t; value: TH
-        weights=model.param_arc_inv_sos1_weights_th,  # key: TH; value: weight
-        sos=1,
+        model.set_T_sos1,
+        rule=rule_constr_arc_group_sos1,
+        sos=1
     )
 
     # *************************************************************************
diff --git a/src/topupopt/problems/esipp/problem.py b/src/topupopt/problems/esipp/problem.py
index ced3fb7209d23097b9ed7161f22c4dd014b5be9a..ce148d327e2d2a21fdd28340d1e02fd6d6d1fd01 100644
--- a/src/topupopt/problems/esipp/problem.py
+++ b/src/topupopt/problems/esipp/problem.py
@@ -3133,12 +3133,6 @@ class InfrastructurePlanningProblem(EnergySystem):
 
         set_GLLJH_arc_inv_sos1 = tuple(param_arc_inv_sos1_weights_glljh.keys())
 
-        set_GLLJH_arc_inv_sos1_gllj = {
-            (g, u, v, j): tuple((g, u, v, j, h) for h in set_H_gllj[(g, u, v, j)])
-            for (g, u, v) in set_J_arc_sos1
-            for j in set_J_arc_sos1[(g, u, v)]
-        }
-
         # sos1 weights for flow sense determination
 
         param_arc_sns_sos1_weights_glljqk = {
@@ -3423,7 +3417,6 @@ class InfrastructurePlanningProblem(EnergySystem):
                 "set_GLLJ_pre_fin_red": set_GLLJ_pre_fin_red,
                 "set_GLLJ_arc_inv_sos1": set_GLLJ_arc_inv_sos1,
                 "set_GLLJH_arc_inv_sos1": set_GLLJH_arc_inv_sos1,
-                "set_GLLJH_arc_inv_sos1_gllj": set_GLLJH_arc_inv_sos1_gllj,
                 "set_GLLJQK_und_sns_sos1_red": set_GLLJQK_und_sns_sos1_red,
                 "set_GLLJ_int": set_GLLJ_int,
                 "set_GLLJ_static": set_GLLJ_static,
diff --git a/tests/test_esipp_problem.py b/tests/test_esipp_problem.py
index 3289150f82a3413f4fda276985de9dacb73d3b6c..580d38d7a87114ce5006bab6dc8046762f4d1622 100644
--- a/tests/test_esipp_problem.py
+++ b/tests/test_esipp_problem.py
@@ -2113,6 +2113,8 @@ class TestESIPPProblem:
     # *************************************************************************
     # *************************************************************************
     
+    # TODO: test arc groups with sos
+    
     def test_arc_groups_individual_ref(self):
         
         # time frame