diff --git a/Project.toml b/Project.toml
index e49f22503..2cd441e3d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -25,7 +25,9 @@ SumTypes = "8e1ec7a9-0e02-4297-b0fe-6433085c89f2"
 [weakdeps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Hecke = "3e1990a7-5d81-5526-99ce-9ba3ff248f21"
+GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 LDPCDecoders = "3c486d74-64b9-4c60-8b1a-13a564e77efb"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PyQDecoders = "17f5de1a-9b79-4409-a58d-4d45812840f7"
@@ -35,6 +37,7 @@ QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678"
 [extensions]
 QuantumCliffordGPUExt = "CUDA"
 QuantumCliffordHeckeExt = "Hecke"
+QuantumCliffordJuMPExt = ["GLPK", "JuMP"]
 QuantumCliffordLDPCDecodersExt = "LDPCDecoders"
 QuantumCliffordMakieExt = "Makie"
 QuantumCliffordPlotsExt = "Plots"
@@ -48,10 +51,12 @@ Combinatorics = "1.0"
 DataStructures = "0.18"
 DocStringExtensions = "0.9"
 Graphs = "1.9"
+GLPK = "1.2"
 Hecke = "0.28, 0.29, 0.30, 0.31, 0.32, 0.33, 0.34"
 HostCPUFeatures = "0.1.6"
 ILog2 = "0.2.3, 1, 2"
 InteractiveUtils = "1.9"
+JuMP = "1.23"
 LDPCDecoders = "0.3.1"
 LinearAlgebra = "1.9"
 MacroTools = "0.5.9"
diff --git a/docs/Project.toml b/docs/Project.toml
index 8ef02232a..ea4ca3714 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -3,10 +3,12 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
+GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 Hecke = "3e1990a7-5d81-5526-99ce-9ba3ff248f21"
 LDPCDecoders = "3c486d74-64b9-4c60-8b1a-13a564e77efb"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PyQDecoders = "17f5de1a-9b79-4409-a58d-4d45812840f7"
diff --git a/docs/src/references.bib b/docs/src/references.bib
index 05276ff75..53902b390 100644
--- a/docs/src/references.bib
+++ b/docs/src/references.bib
@@ -561,3 +561,10 @@ @article{haah2011local
   pages={042330},
   year={2011},
 }
+
+@article{landahl2011fault,
+  title={Fault-tolerant quantum computing with color codes},
+  author={Landahl, Andrew J and Anderson, Jonas T and Rice, Patrick R},
+  journal={arXiv preprint arXiv:1108.5738},
+  year={2011}
+}
diff --git a/ext/QuantumCliffordHeckeExt/QuantumCliffordHeckeExt.jl b/ext/QuantumCliffordHeckeExt/QuantumCliffordHeckeExt.jl
index 89d7bae0c..5fb858826 100644
--- a/ext/QuantumCliffordHeckeExt/QuantumCliffordHeckeExt.jl
+++ b/ext/QuantumCliffordHeckeExt/QuantumCliffordHeckeExt.jl
@@ -12,7 +12,7 @@ import Nemo: characteristic, matrix_repr, GF, ZZ, lift
 import QuantumClifford.ECC: AbstractECC, CSS, ClassicalCode,
     hgp, code_k, code_n, code_s, iscss, parity_checks, parity_checks_x, parity_checks_z, parity_checks_xz,
     two_block_group_algebra_codes, generalized_bicycle_codes, bicycle_codes, check_repr_commutation_relation,
-    haah_cubic_codes
+    haah_cubic_codes, minimum_distance
 
 include("util.jl")
 include("types.jl")
diff --git a/ext/QuantumCliffordHeckeExt/lifted_product.jl b/ext/QuantumCliffordHeckeExt/lifted_product.jl
index 66bd35ed6..9f17c5963 100644
--- a/ext/QuantumCliffordHeckeExt/lifted_product.jl
+++ b/ext/QuantumCliffordHeckeExt/lifted_product.jl
@@ -162,7 +162,9 @@ Here is an example of a [[56, 28, 2]] 2BGA code from Table 2 of [lin2024quantum]
 with direct product of `C₄ x C₂`.
 
 ```jldoctest
-julia> import Hecke: group_algebra, GF, abelian_group, gens
+julia> import Hecke: group_algebra, GF, abelian_group, gens; import GLPK; import JuMP;
+
+julia> using QuantumClifford.ECC: minimum_distance, two_block_group_algebra_codes; # hide
 
 julia> GA = group_algebra(GF(2), abelian_group([14,2]));
 
@@ -174,8 +176,8 @@ julia> B = 1 + x^7 + s + x^8 + s*x^7 + x
 
 julia> c = two_block_group_algebra_codes(A,B);
 
-julia> code_n(c), code_k(c)
-(56, 28)
+julia> code_n(c), code_k(c), minimum_distance(c)
+(56, 28, 2)
 ```
 
 ### Bivariate Bicycle codes
@@ -189,7 +191,7 @@ The ECC Zoo has an [entry for this family](https://errorcorrectionzoo.org/c/qcga
 A [[756, 16, ≤ 34]] code from Table 3 of [bravyi2024high](@cite):
 
 ```jldoctest
-julia> import Hecke: group_algebra, GF, abelian_group, gens
+julia> import Hecke: group_algebra, GF, abelian_group, gens;
 
 julia> l=21; m=18;
 
@@ -215,10 +217,16 @@ where `𝐺ᵣ = ℤ/l₁ × ℤ/l₂ × ... × ℤ/lᵣ`.
 A [[48, 4, 6]] Weight-6 TB-QLDPC code from Appendix A Table 2 of [voss2024multivariatebicyclecodes](@cite).
 
 ```jldoctest
-julia> import Hecke: group_algebra, GF, abelian_group, gens; # hide
+julia> import Hecke: group_algebra, GF, abelian_group, gens; import GLPK; import JuMP;
+
+julia> using QuantumClifford.ECC: minimum_distance, two_block_group_algebra_codes; # hide
 
 julia> l=4; m=6;
 
+julia> GA = group_algebra(GF(2), abelian_group([l, m]));
+
+julia> x, y = gens(GA);
+
 julia> z = x*y;
 
 julia> A = x^3 + y^5;
@@ -227,8 +235,8 @@ julia> B = x + z^5 + y^5 + y^2;
 
 julia> c = two_block_group_algebra_codes(A, B);
 
-julia> code_n(c), code_k(c)
-(48, 4)
+julia> code_n(c), code_k(c), minimum_distance(c)
+(48, 4, 6)
 ```
 
 ### Coprime Bivariate Bicycle code
@@ -241,7 +249,9 @@ based on abelian group `ℤₗ x ℤₘ` where `ℤⱼ` cyclic group of order `j
 [108, 12, 6]] coprime-bivariate bicycle (BB) code from Table 2 of [wang2024coprime](@cite).
 
 ```jldoctest
-julia> import Hecke: group_algebra, GF, abelian_group, gens;
+julia> import Hecke: group_algebra, GF, abelian_group, gens; import GLPK; import JuMP;
+
+julia> using QuantumClifford.ECC: minimum_distance, two_block_group_algebra_codes; # hide
 
 julia> l=2; m=27;
 
@@ -252,7 +262,11 @@ julia> 𝜋 = gens(GA)[1];
 julia> A = 𝜋^2 + 𝜋^5  + 𝜋^44;
 
 julia> B = 𝜋^8 + 𝜋^14 + 𝜋^47;
-(108, 12)
+
+julia> c = two_block_group_algebra_codes(A, B);
+
+julia> code_n(c), code_k(c), minimum_distance(c)
+(108, 12, 6)
 ```
 
 See also: [`LPCode`](@ref), [`generalized_bicycle_codes`](@ref), [`bicycle_codes`](@ref), [`haah_cubic_codes`](@ref).
diff --git a/ext/QuantumCliffordJuMPExt/QuantumCliffordJuMPExt.jl b/ext/QuantumCliffordJuMPExt/QuantumCliffordJuMPExt.jl
new file mode 100644
index 000000000..512137318
--- /dev/null
+++ b/ext/QuantumCliffordJuMPExt/QuantumCliffordJuMPExt.jl
@@ -0,0 +1,24 @@
+module QuantumCliffordJuMPExt
+
+using DocStringExtensions
+
+import JuMP
+import JuMP: @variable, @objective, @constraint, optimize!, Model, set_silent,
+    sum, value
+import GLPK
+import GLPK: Optimizer
+
+import QuantumClifford
+import QuantumClifford: stab_to_gf2, logicalxview, logicalzview, canonicalize!,
+    MixedDestabilizer, Stabilizer
+import QuantumClifford.ECC
+import QuantumClifford.ECC: minimum_distance, AbstractECC, code_n, code_k,
+    parity_checks
+
+import SparseArrays
+import SparseArrays: SparseMatrixCSC, sparse
+
+include("util.jl")
+include("min_distance_mixed_integer_programming.jl")
+
+end
diff --git a/ext/QuantumCliffordJuMPExt/min_distance_mixed_integer_programming.jl b/ext/QuantumCliffordJuMPExt/min_distance_mixed_integer_programming.jl
new file mode 100644
index 000000000..fd6bd5cb0
--- /dev/null
+++ b/ext/QuantumCliffordJuMPExt/min_distance_mixed_integer_programming.jl
@@ -0,0 +1,147 @@
+"""
+TYPEDSIGNATURES
+
+Computes the minimum Hamming weight of a binary vector `x` by solving an **mixed
+integer program (MIP)** that satisfies the following constraints:
+
+- \$\\text{stab} \\cdot x \\equiv 0 \\pmod{2}\$: The binary vector x must have an
+even overlap with each X-check of the stabilizer binary representation `stab`.
+- \$\\text{logicOp} \\cdot x \\equiv 1 \\pmod{2}\$: The binary vector x must have
+an odd overlap with logical-X operator `logicOp` on the i-th logical qubit.
+
+It computes the minimum Hamming weight \$d_{Z}\$ for the Z-type logical
+operator. The distance for X-type logical operators is the same.
+
+### Mixed Integer Programming
+
+The MIP minimizes the Hamming weight `w(x)`, defined as the number of nonzero
+elements in `x`, while satisfying the constraints:
+
+\$\$
+\\begin{aligned}
+    \\text{Minimize} \\quad & w(x) = \\sum_{i=1}^n x_i, \\\\
+    \\text{subject to} \\quad & \\text{stab} \\cdot x \\equiv 0 \\pmod{2}, \\\\
+                           & \\text{logicOp} \\cdot x \\equiv 1 \\pmod{2}, \\\\
+                           & x_i \\in \\{0, 1\\} \\quad \\text{for all } i.
+\\end{aligned}
+\$\$
+
+Here:
+- \$\\text{stab}\$ is the binary matrix representing the stabilizer group.
+- \$\\text{logicOp}\$ is the binary vector representing the logical-X operator.
+- `x` is the binary vector (decision variable) being optimized.
+
+The optimal solution \$w_{i}\$ for each logical-X operator corresponds to the
+minimum weight of a Pauli Z-type operator satisfying the above conditions.
+The Z-type distance is given by:
+
+\$\$
+\\begin{aligned}
+    d_Z = \\min(w_1, w_2, \\dots, w_k),
+\\end{aligned}
+\$\$
+
+where k is the number of logical qubits.
+
+# Example
+
+A [[40, 8, 5]] 2BGA code with the minimum distance of 5.
+
+```jldoctest
+julia> import Hecke: group_algebra, GF, abelian_group, gens; import GLPK; import JuMP;
+
+julia> using QuantumClifford.ECC: minimum_distance, two_block_group_algebra_codes;
+
+julia> l = 10; m = 2;
+
+julia> GA = group_algebra(GF(2), abelian_group([l,m]));
+
+julia> x, s = gens(GA);
+
+julia> A = 1 + x^6;
+
+julia> B = 1 + x^5 + s + x^6 + x + s*x^2;
+
+julia> c = two_block_group_algebra_codes(A,B);
+
+julia> minimum_distance(c)
+5
+```
+
+### Applications
+
+- The first usecase of the MIP approach was the code capacity Most Likely
+Error (MLE) decoder for color codes introduced in [landahl2011fault](@cite).
+- For all quantum LDPC codes presented in [panteleev2021degenerate](@cite),
+the lower and upper bounds on the minimum distance was obtained by reduction
+to a mixed integer linear program and using the GNU Linear Programming Kit.
+- For all the Bivariate Bicycle (BB) codes presented in [bravyi2024high](@cite),
+the code distance was calculated using the mixed integer programming approach.
+
+"""
+function minimum_distance(c::Stabilizer; num_logical_qubits=code_k(c))
+    lx, _ = get_lx_lz(c)
+    H = stab_to_gf2(parity_checks(c))
+    H = sparse(H)
+    hx =  Matrix{Int}(get_stab_hx(H))
+    num_logical_qubits == 1 && return _minimum_distance(hx, lx[num_logical_qubits, :])
+    if 2 <= num_logical_qubits && num_logical_qubits <= code_k(c)
+        w_values = []
+        for i in 1:num_logical_qubits
+            w = _minimum_distance(hx, lx[num_logical_qubits, :])
+            push!(w_values, w)
+        end
+        return round(Int, sum(w_values)/num_logical_qubits)
+    else
+        throw(ArgumentError("The number of logical qubits must be between 1 to $(code_k(c)) inclusive"))
+    end
+end
+
+# Computing minimum distance for quantum LDPC codes is a NP-Hard problem,
+# but we can solve mixed integer program (MIP) for small instance sizes.
+# inspired from [bravyi2024high](@cite) and [landahl2011fault](@cite).
+function _minimum_distance(hx, lx)
+    n = size(hx, 2) # number of qubits
+    m = size(hx, 1) # number of stabilizers
+    # Maximum stabilizer weight
+    whx = maximum(sum(hx[i, :] for i in 1:m))  # Maximum row sum of hx
+    # Weight of the logical operator
+    wlx = count(!iszero, lx)  # Count non-zero elements in lx
+    # Number of slack variables needed to express orthogonality constraints modulo two
+    num_anc_hx = ceil(Int, log2(whx))
+    num_anc_logical = ceil(Int, log2(wlx))
+    num_var = n + m * num_anc_hx + num_anc_logical
+    model = Model(GLPK.Optimizer) # model
+    set_silent(model)
+    @variable(model, x[1:num_var], Bin) # binary variables
+    @objective(model, Min, sum(x[i] for i in 1:n)) # objective function
+    # Orthogonality to rows of hx constraints
+    for row in 1:m
+        weight = zeros(Int, num_var)
+        supp = findall(hx[row, :] .!= 0)  # Non-zero indices in hx[row, :]
+             for q in supp
+                 weight[q] = 1
+             end
+             cnt = 1
+             for q in 1:num_anc_hx
+                 weight[n + (row - 1) * num_anc_hx + q] = -(1 << cnt)
+             cnt += 1
+         end
+        @constraint(model, sum(weight[i] * x[i] for i in 1:num_var) == 0)
+    end
+    # Odd overlap with lx constraint
+    supp = findall(lx .!= 0)  # Non-zero indices in lx
+    weight = zeros(Int, num_var)
+    for q in supp
+        weight[q] = 1
+    end
+    cnt = 1
+    for q in 1:num_anc_logical
+        weight[n + m * num_anc_hx + q] = -(1 << cnt)
+        cnt += 1
+    end
+    @constraint(model, sum(weight[i] * x[i] for i in 1:num_var) == 1)
+    optimize!(model)
+    opt_val = sum(value(x[i]) for i in 1:n)
+    return Int(opt_val)
+end
diff --git a/ext/QuantumCliffordJuMPExt/util.jl b/ext/QuantumCliffordJuMPExt/util.jl
new file mode 100644
index 000000000..8626f4695
--- /dev/null
+++ b/ext/QuantumCliffordJuMPExt/util.jl
@@ -0,0 +1,14 @@
+function get_stab_hx(matrix::SparseMatrixCSC{T, Int}) where T
+    rows, cols = size(matrix)
+    rhs_start = div(cols, 2) + 1
+    rhs_cols = matrix[:, rhs_start:cols]
+    non_zero_rows_rhs = unique(rhs_cols.rowval)
+    zero_rows_rhs = setdiff(1:rows, non_zero_rows_rhs)
+    return matrix[zero_rows_rhs, :]
+end
+
+function get_lx_lz(c::Stabilizer)
+    lx = stab_to_gf2(logicalxview(canonicalize!(MixedDestabilizer(c))))
+    lz = stab_to_gf2(logicalzview(canonicalize!(MixedDestabilizer(c))))
+    return lx, lz
+end
diff --git a/src/ecc/ECC.jl b/src/ecc/ECC.jl
index 231816318..8f36fa06a 100644
--- a/src/ecc/ECC.jl
+++ b/src/ecc/ECC.jl
@@ -21,7 +21,7 @@ using Nemo: ZZ, residue_ring, matrix, finite_field, GF, minpoly, coeff, lcm, FqP
 abstract type AbstractECC end
 
 export parity_checks, parity_checks_x, parity_checks_z, iscss,
-    code_n, code_s, code_k, rate, distance,
+    code_n, code_s, code_k, rate, distance, minimum_distance,
     isdegenerate, faults_matrix,
     naive_syndrome_circuit, shor_syndrome_circuit, naive_encoding_circuit,
     RepCode, LiftedCode,
@@ -123,6 +123,16 @@ end
 """The distance of a code."""
 function distance end
 
+"""The minimum distance of a qLDPC code.
+
+Requires a JuMP.jl backend to be loaded, e.g. `using JuMP`.
+
+Requires a GLPK.jl backend to be loaded, e.g. `using GLPK`.
+"""
+function minimum_distance end
+
+minimum_distance(c::AbstractECC; kwargs...) = minimum_distance(parity_checks(c); kwargs...)
+
 """Parity matrix of a code, given as a stabilizer tableau."""
 function parity_matrix(c::AbstractECC)
     paritym = stab_to_gf2(parity_checks(c::AbstractECC))