diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml
index bba03cd75..4fa9f2804 100644
--- a/.github/workflows/ci.yaml
+++ b/.github/workflows/ci.yaml
@@ -53,16 +53,6 @@ jobs:
       run: |
         echo -ne "url: ${CDSAPI_URL}\nkey: ${CDSAPI_TOKEN}\n" > ~/.cdsapirc
 
-    - name: Add solver to environment
-      run: |
-        echo -e "- glpk\n- ipopt<3.13.3" >> envs/environment.yaml
-      if: ${{ matrix.os }} == 'windows-latest'
-
-    - name: Add solver to environment
-      run: |
-        echo -e "- glpk\n- ipopt" >> envs/environment.yaml
-      if: ${{ matrix.os }} != 'windows-latest'
-
     - name: Setup micromamba
       uses: mamba-org/setup-micromamba@v1
       with:
diff --git a/doc/release_notes.rst b/doc/release_notes.rst
index 867632874..6fa887385 100644
--- a/doc/release_notes.rst
+++ b/doc/release_notes.rst
@@ -60,6 +60,11 @@ Upcoming Release
 * The rule ``plot_network`` has been split into separate rules for plotting
   electricity, hydrogen and gas networks.
 
+* To determine the optimal topology to meet the number of clusters, the workflow used pyomo in combination with ``ipopt`` or ``gurobi``. This dependency has been replaced by using ``linopy`` in combination with ``scipopt`` or ``gurobi``. The environment file has been updated accordingly.
+
+* The ``highs`` solver was added to the default environment file.
+
+
 
 PyPSA-Eur 0.9.0 (5th January 2024)
 ==================================
diff --git a/envs/environment.yaml b/envs/environment.yaml
index 6ff4b7f15..26e18f0d0 100644
--- a/envs/environment.yaml
+++ b/envs/environment.yaml
@@ -35,8 +35,9 @@ dependencies:
 - netcdf4
 - networkx
 - scipy
+- glpk
 - shapely>=2.0
-- pyomo
+- pyscipopt
 - matplotlib
 - proj
 - fiona
@@ -47,7 +48,6 @@ dependencies:
 - tabula-py
 - pyxlsb
 - graphviz
-- ipopt
 
   # Keep in conda environment when calling ipython
 - ipython
@@ -60,3 +60,4 @@ dependencies:
 
 - pip:
   - tsam>=2.3.1
+  - highspy
diff --git a/scripts/cluster_network.py b/scripts/cluster_network.py
index 01af29aa5..39b6ad966 100644
--- a/scripts/cluster_network.py
+++ b/scripts/cluster_network.py
@@ -122,14 +122,15 @@
 """
 
 import logging
+import os
 import warnings
 from functools import reduce
 
 import geopandas as gpd
+import linopy
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
-import pyomo.environ as po
 import pypsa
 import seaborn as sns
 from _helpers import configure_logging, update_p_nom_max
@@ -214,7 +215,7 @@ def get_feature_for_hac(n, buses_i=None, feature=None):
     return feature_data
 
 
-def distribute_clusters(n, n_clusters, focus_weights=None, solver_name="cbc"):
+def distribute_clusters(n, n_clusters, focus_weights=None, solver_name="scip"):
     """
     Determine the number of clusters per country.
     """
@@ -254,31 +255,22 @@ def distribute_clusters(n, n_clusters, focus_weights=None, solver_name="cbc"):
         L.sum(), 1.0, rtol=1e-3
     ), f"Country weights L must sum up to 1.0 when distributing clusters. Is {L.sum()}."
 
-    m = po.ConcreteModel()
-
-    def n_bounds(model, *n_id):
-        return (1, N[n_id])
-
-    m.n = po.Var(list(L.index), bounds=n_bounds, domain=po.Integers)
-    m.tot = po.Constraint(expr=(po.summation(m.n) == n_clusters))
-    m.objective = po.Objective(
-        expr=sum((m.n[i] - L.loc[i] * n_clusters) ** 2 for i in L.index),
-        sense=po.minimize,
+    m = linopy.Model()
+    clusters = m.add_variables(
+        lower=1, upper=N, coords=[L.index], name="n", integer=True
     )
-
-    opt = po.SolverFactory(solver_name)
-    if solver_name == "appsi_highs" or not opt.has_capability("quadratic_objective"):
-        logger.warning(
-            f"The configured solver `{solver_name}` does not support quadratic objectives. Falling back to `ipopt`."
+    m.add_constraints(clusters.sum() == n_clusters, name="tot")
+    # leave out constant in objective (L * n_clusters) ** 2
+    m.objective = (clusters * clusters - 2 * clusters * L * n_clusters).sum()
+    if solver_name == "gurobi":
+        logging.getLogger("gurobipy").propagate = False
+    elif solver_name != "scip":
+        logger.info(
+            f"The configured solver `{solver_name}` does not support quadratic objectives. Falling back to `scip`."
         )
-        opt = po.SolverFactory("ipopt")
-
-    results = opt.solve(m)
-    assert (
-        results["Solver"][0]["Status"] == "ok"
-    ), f"Solver returned non-optimally: {results}"
-
-    return pd.Series(m.n.get_values(), index=L.index).round().astype(int)
+        solver_name = "scip"
+    m.solve(solver_name=solver_name)
+    return m.solution["n"].to_series().astype(int)
 
 
 def busmap_for_n_clusters(
@@ -372,7 +364,7 @@ def busmap_for_country(x):
 
     return (
         n.buses.groupby(["country", "sub_network"], group_keys=False)
-        .apply(busmap_for_country)
+        .apply(busmap_for_country, include_groups=False)
         .squeeze()
         .rename("busmap")
     )
@@ -385,7 +377,7 @@ def clustering_for_n_clusters(
     aggregate_carriers=None,
     line_length_factor=1.25,
     aggregation_strategies=dict(),
-    solver_name="cbc",
+    solver_name="scip",
     algorithm="hac",
     feature=None,
     extended_link_costs=0,
@@ -462,7 +454,6 @@ def plot_busmap_for_n_clusters(n, n_clusters, fn=None):
 
     params = snakemake.params
     solver_name = snakemake.config["solving"]["solver"]["name"]
-    solver_name = "appsi_highs" if solver_name == "highs" else solver_name
 
     n = pypsa.Network(snakemake.input.network)