diff --git a/src/utilities.jl b/src/utilities.jl
index c3aa6f1..6c9b8e5 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -6,4 +6,272 @@ function _copy_model(
     model::M
     ) where {M <: JuMP.AbstractModel}
     return M()
-end
\ No newline at end of file
+end
+
+"""
+    JuMP.copy_extension_data(
+        data::GDPData,
+        new_model::JuMP.AbstractModel,
+        old_model::JuMP.AbstractModel
+    )::GDPData
+
+Extend `JuMP.copy_extension_data` to initialize an empty [`GDPData`](@ref) object 
+for the copied model. This is the first step in the model copying process and is 
+automatically called by `JuMP.copy_model`. The actual GDP data (logical variables, 
+disjunctions, etc.) is copied separately via [`copy_gdp_data`](@ref).
+"""
+function JuMP.copy_extension_data(
+    data::GDPData{M, V, C, T},
+    new_model::JuMP.AbstractModel,
+    old_model::JuMP.AbstractModel
+) where {M, V, C, T}
+    return GDPData{M, V, C}()
+end
+
+"""
+    copy_gdp_data(
+        model::JuMP.AbstractModel,
+        new_model::JuMP.AbstractModel,
+        ref_map::JuMP.GenericReferenceMap
+    )::Dict{LogicalVariableRef, LogicalVariableRef}
+
+Copy all GDP-specific data from `model` to `new_model`, including logical variables, 
+logical constraints, disjunct constraints, and disjunctions. This function is called 
+automatically by [`copy_gdp_model`](@ref) after `JuMP.copy_model` has copied the base 
+model structure.
+
+**Arguments**
+- `model::JuMP.AbstractModel`: The source model containing GDP data to copy.
+- `new_model::JuMP.AbstractModel`: The destination model that will receive the copied GDP data.
+- `ref_map::JuMP.GenericReferenceMap`: The reference map from `JuMP.copy_model` that maps 
+  old variable references to new ones.
+
+**Returns**
+- `Dict{LogicalVariableRef, LogicalVariableRef}`: A mapping from old logical variable 
+  references to new logical variable references.
+"""
+function copy_gdp_data(
+    model::M,
+    new_model::M,
+    ref_map::GenericReferenceMap
+    ) where {M <: JuMP.AbstractModel}
+    
+    old_gdp = model.ext[:GDP]
+
+    # GDPData contains the following fields.
+    # DICTIONARIES (for loops below)
+    # - logical_variables
+    # - logical_constraints
+    # - disjunct_constraints
+    # - disjunctions
+    # - exactly1_constraints
+    # - indicator_to_binary
+    # - indicator_to_constraints
+    # - constraint_to_indicator
+    # - variable_bounds
+    # SINGLE VALUES (copy directly)
+    # - solution_method
+    # - ready_to_optimize
+
+    new_gdp = new_model.ext[:GDP]
+
+    # Creating maps from old to new model.
+    var_map = Dict(v => ref_map[v] for v in all_variables(model))
+    lv_map = Dict{LogicalVariableRef{M}, LogicalVariableRef{M}}()
+    lc_map = Dict{LogicalConstraintRef{M}, LogicalConstraintRef{M}}()
+    disj_map = Dict{DisjunctionRef{M}, DisjunctionRef{M}}()
+    disj_con_map = Dict{DisjunctConstraintRef{M}, DisjunctConstraintRef{M}}()
+
+    # Copying logical variables
+    for (idx, var) in old_gdp.logical_variables
+        old_var_ref = LogicalVariableRef(model, idx)
+        new_var_data = LogicalVariableData(var.variable, var.name)
+        new_var = LogicalVariableRef(new_model, idx)
+        lv_map[old_var_ref] = new_var
+        # Update to new_gdp.logical_variables
+        new_gdp.logical_variables[idx] = new_var_data
+    end
+
+    # Copying logical constraints
+    for (idx, lc_data) in old_gdp.logical_constraints
+        old_con_ref = LogicalConstraintRef(model, idx)
+        new_con_ref = LogicalConstraintRef(new_model, idx)
+        c = lc_data.constraint
+        expr = _replace_variables_in_constraint(c.func, lv_map)
+        new_con = JuMP.build_constraint(error, expr, c.set)
+        JuMP.add_constraint(new_model, new_con, lc_data.name)
+        lc_map[old_con_ref] = new_con_ref
+    end
+
+    # Copying disjunct constraints
+    for (idx, disj_con_data) in old_gdp.disjunct_constraints
+        old_constraint = disj_con_data.constraint
+        old_dc_ref = DisjunctConstraintRef(model, idx)
+        old_indicator = old_gdp.constraint_to_indicator[old_dc_ref]
+        new_indicator = lv_map[old_indicator]
+        new_expr = _replace_variables_in_constraint(old_constraint.func, 
+        var_map
+        )
+        # Update to new_gdp.disjunct_constraints
+        new_con = JuMP.build_constraint(error, new_expr, 
+        old_constraint.set, Disjunct(new_indicator)
+        )
+        new_dc_ref = JuMP.add_constraint(new_model, new_con, disj_con_data.name)
+        disj_con_map[old_dc_ref] = new_dc_ref
+    end
+
+    # Copying disjunctions
+    for (idx, disj_data) in old_gdp.disjunctions
+        old_disj = disj_data.constraint
+        new_indicators = [_replace_variables_in_constraint(indicator, lv_map) 
+        for indicator in old_disj.indicators
+            ]
+        new_disj = Disjunction(new_indicators, old_disj.nested)
+        disj_map[DisjunctionRef(model, idx)] = DisjunctionRef(new_model, idx)
+        # Update to new_gdp.disjunctions
+        new_gdp.disjunctions[idx] = ConstraintData(new_disj, disj_data.name)
+    end
+
+    # Copying exactly1 constraints
+    for (d_ref, lc_ref) in old_gdp.exactly1_constraints
+        new_lc_ref = lc_map[lc_ref]
+        new_d_ref = disj_map[d_ref]
+        # Update to new_gdp.exactly1_constraints
+        new_gdp.exactly1_constraints[new_d_ref] = new_lc_ref
+    end
+
+    # Copying indicator to binary
+    for (lv_ref, bref) in old_gdp.indicator_to_binary
+        new_bref = _remap_indicator_to_binary(bref, var_map)
+        # Update to new_gdp.indicator_to_binary
+        new_gdp.indicator_to_binary[lv_map[lv_ref]] = new_bref
+    end
+
+    # Copying indicator to constraints
+    for (lv_ref, con_refs) in old_gdp.indicator_to_constraints
+        new_lvar_ref = lv_map[lv_ref]
+        new_con_refs = Vector{Union{DisjunctConstraintRef{M}, DisjunctionRef{M}}}()
+        for con_ref in con_refs
+            new_con_ref = _remap_indicator_to_constraint(con_ref, 
+            disj_con_map, disj_map
+            )
+            push!(new_con_refs, new_con_ref)
+        end
+        # Update to new_gdp.indicator_to_constraints
+        new_gdp.indicator_to_constraints[new_lvar_ref] = new_con_refs
+    end
+
+    # Copying constraint to indicator
+    for (con_ref, lv_ref) in old_gdp.constraint_to_indicator
+        # Update to new_gdp.constraint_to_indicator
+        new_gdp.constraint_to_indicator[
+            _remap_constraint_to_indicator(con_ref, disj_con_map, disj_map)
+            ] = lv_map[lv_ref]
+    end
+
+    # Copying variable bounds
+    for (v, bounds) in old_gdp.variable_bounds
+        # Update to new_gdp.variable_bounds
+        new_gdp.variable_bounds[var_map[v]] = bounds
+    end
+
+    # Copying solution method and ready to optimize
+    new_gdp.solution_method = old_gdp.solution_method
+    new_gdp.ready_to_optimize = old_gdp.ready_to_optimize
+
+    return lv_map
+end
+
+"""
+    copy_gdp_model(model::JuMP.AbstractModel)
+
+Create a copy of a [`GDPModel`](@ref), including all variables, constraints, and 
+GDP-specific data (logical variables, disjunctions, etc.).
+
+**Arguments**
+- `model::JuMP.AbstractModel`: The GDP model to copy.
+
+**Returns**
+A tuple `(new_model, ref_map, lv_map)` where:
+- `new_model`: The copied model.
+- `ref_map::JuMP.GenericReferenceMap`: Maps old variable and constraint references to new ones.
+- `lv_map::Dict{LogicalVariableRef, LogicalVariableRef}`: Maps old logical variable 
+  references to new ones.
+
+## Example
+```julia
+using DisjunctiveProgramming, HiGHS
+model = GDPModel(HiGHS.Optimizer)
+@variable(model, x)
+@variable(model, Y[1:2], LogicalVariable)
+@constraint(model, x <= 10, Disjunct(Y[1]))
+@constraint(model, x >= 20, Disjunct(Y[2]))
+@disjunction(model, Y)
+
+new_model, ref_map, lv_map = copy_gdp_model(model)
+```
+"""
+function copy_gdp_model(model::M) where {M <: JuMP.AbstractModel}
+    new_model, ref_map = JuMP.copy_model(model)
+    lv_map = copy_gdp_data(model, new_model, ref_map)
+    return new_model, ref_map, lv_map
+end
+################################################################################
+#                                GDP REMAPPING
+################################################################################
+# These remapping functions use multiple dispatch to handle different types that
+# can appear in GDP data structures during model copying.
+#
+# Indicators can be represented by a variable or an affine expression to 
+# indicate a complementary relationship with another variable.
+# This translates to a binary or affine expression in its binary reformulation.
+#
+# Depending on the above, different mappings are required for indicator_to_binary,
+# indicator_to_constraints, and constraint_to_indicator.
+################################################################################
+
+function _remap_indicator_to_constraint(
+    con_ref::DisjunctConstraintRef,
+    disj_con_map::Dict{DisjunctConstraintRef{M}, DisjunctConstraintRef{M}},
+    ::Dict{DisjunctionRef{M}, DisjunctionRef{M}}
+) where {M <: JuMP.AbstractModel}
+    return disj_con_map[con_ref]   
+end
+
+function _remap_indicator_to_constraint(
+    con_ref::DisjunctionRef,
+    ::Dict{DisjunctConstraintRef{M}, DisjunctConstraintRef{M}},
+    disj_map::Dict{DisjunctionRef{M}, DisjunctionRef{M}}
+) where {M <: JuMP.AbstractModel}
+    return disj_map[con_ref]   
+end
+
+function _remap_indicator_to_binary(
+    bref::JuMP.AbstractVariableRef,
+    var_map::Dict{V, V}
+) where {V <: JuMP.AbstractVariableRef}
+    return var_map[bref]
+end
+
+function _remap_indicator_to_binary(
+    bref::JuMP.GenericAffExpr,
+    var_map::Dict{V, V}
+) where {V <: JuMP.AbstractVariableRef}
+    return _replace_variables_in_constraint(bref, var_map)
+end
+
+function _remap_constraint_to_indicator(
+    con_ref::DisjunctConstraintRef,
+    disj_con_map::Dict{DisjunctConstraintRef{M}, DisjunctConstraintRef{M}},
+    ::Dict{DisjunctionRef{M}, DisjunctionRef{M}}
+) where {M <: JuMP.AbstractModel}
+    return disj_con_map[con_ref]   
+end
+
+function _remap_constraint_to_indicator(
+    con_ref::DisjunctionRef,
+    ::Dict{DisjunctConstraintRef{M}, DisjunctConstraintRef{M}},
+    disj_map::Dict{DisjunctionRef{M}, DisjunctionRef{M}}
+) where {M <: JuMP.AbstractModel}
+    return disj_map[con_ref]   
+end
diff --git a/test/model.jl b/test/model.jl
index 1cd2bff..21e7709 100644
--- a/test/model.jl
+++ b/test/model.jl
@@ -41,10 +41,158 @@ function test_set_optimizer()
     @test solver_name(model) == "HiGHS"
 end
 
+function test_remapping_functions()
+    model = GDPModel()
+    @variable(model, x)
+    @variable(model, y)
+    @variable(model, z)
+    
+    var_map = Dict{VariableRef, VariableRef}(x => y, z => x)
+    
+    @test DP._remap_indicator_to_binary(x, var_map) == y
+    @test DP._remap_indicator_to_binary(z, var_map) == x
+    
+    aff_expr = 2.0 * x + 3.0 * z + 5.0
+    remapped_expr = DP._remap_indicator_to_binary(aff_expr, var_map)
+    @test remapped_expr isa JuMP.GenericAffExpr
+    @test remapped_expr.constant == 5.0
+    @test remapped_expr.terms[y] == 2.0  # x remapped to y
+    @test remapped_expr.terms[x] == 3.0  # z remapped to x
+    @test !haskey(remapped_expr.terms, z)  # z shouldn't exist anymore
+    
+    model2 = GDPModel()
+    @variable(model2, a[1:3])
+    @variable(model2, Y[1:3], DP.Logical)
+    
+    con1 = @constraint(model2, a[1] ≤ 5, DP.Disjunct(Y[1]))
+    con2 = @constraint(model2, a[2] ≥ 2, DP.Disjunct(Y[2]))
+    con3 = @constraint(model2, a[3] == 4, DP.Disjunct(Y[3]))
+    
+    disj1 = DP.@disjunction(model2, Y[1:2])
+    disj2 = DP.@disjunction(model2, Y[2:3])
+    
+    disj_con_map = Dict{DisjunctConstraintRef{Model}, DisjunctConstraintRef{Model}}(
+        con1 => con2,
+        con2 => con3
+    )
+    disj_map = Dict{DisjunctionRef{Model}, DisjunctionRef{Model}}()
+    
+    @test DP._remap_constraint_to_indicator(con1, disj_con_map, disj_map) == con2
+    @test DP._remap_constraint_to_indicator(con2, disj_con_map, disj_map) == con3
+    
+    empty_disj_con_map = Dict{DisjunctConstraintRef{Model}, DisjunctConstraintRef{Model}}()
+    disj_map2 = Dict{DisjunctionRef{Model}, DisjunctionRef{Model}}(
+        disj1 => disj2
+    )
+    
+    @test DP._remap_constraint_to_indicator(disj1, empty_disj_con_map, disj_map2) == disj2
+end
+
+function test_copy_model()
+    model = DP.GDPModel(HiGHS.Optimizer)
+    @variable(model, 0 ≤ x[1:2] ≤ 20)
+    @variable(model, Y[1:2], DP.Logical)
+    @constraint(model, [i = 1:2], [2,5][i] ≤ x[i] ≤ [6,9][i], DP.Disjunct(Y[1]))
+    @constraint(model, [i = 1:2], [8,10][i] ≤ x[i] ≤ [11,15][i], DP.Disjunct(Y[2]))
+    DP.@disjunction(model, Y)
+    DP._variable_bounds(model)[x[1]] = set_variable_bound_info(x[1], BigM())
+    DP._variable_bounds(model)[x[2]] = set_variable_bound_info(x[2], BigM())
+    
+    new_model, ref_map = JuMP.copy_model(model)
+    @test haskey(new_model.ext, :GDP)
+    lv_map = DP.copy_gdp_data(model, new_model, ref_map)
+    @test length(lv_map) == 2 
+    new_model1, ref_map1, lv_map1 = DP.copy_gdp_model(model)
+    @test haskey(new_model1.ext, :GDP)
+    @test length(lv_map1) == 2
+
+    orig_num_vars = num_variables(model)
+    orig_num_constrs = num_constraints(model; 
+    count_variable_in_set_constraints = false
+    )
+    orig1_num_vars = num_variables(new_model1)
+    orig1_num_constrs = num_constraints(new_model1; 
+    count_variable_in_set_constraints = false
+    )
+    
+    @variable(new_model, z >= 0)
+    @constraint(new_model, z <= 100)
+    
+    @test num_variables(model) == orig_num_vars
+    num_con_m = num_constraints(model; 
+    count_variable_in_set_constraints = false)
+    num_con_m1 = num_constraints(new_model1; 
+    count_variable_in_set_constraints = false)
+    @test num_con_m == orig_num_constrs
+    @test num_con_m1 == orig1_num_constrs
+    @test !haskey(object_dictionary(model), :z)
+    
+    @test num_variables(new_model1) == orig1_num_vars
+    
+    @test !haskey(object_dictionary(new_model1), :z)
+
+    @test num_variables(new_model) == orig_num_vars + 1
+    num_con_m2 = num_constraints(new_model; 
+    count_variable_in_set_constraints = false)
+    @test num_con_m2 == orig_num_constrs + 1
+    
+
+    orig_num_lvars = length(DP._logical_variables(model))
+    orig_num_disj_cons = length(DP._disjunct_constraints(model))
+    orig_num_disj = length(DP._disjunctions(model))
+    orig1_num_lvars = length(DP._logical_variables(new_model1))
+    orig1_num_disj_cons = length(DP._disjunct_constraints(new_model1))
+    orig1_num_disj = length(DP._disjunctions(new_model1))
+    
+    @variable(new_model, W[1:2], DP.Logical)
+    @constraint(new_model, z >= 5, DP.Disjunct(W[1]))
+    @constraint(new_model, z <= 3, DP.Disjunct(W[2]))
+    DP.@disjunction(new_model, W)
+    
+    @test length(DP._logical_variables(model)) == orig_num_lvars
+    @test length(DP._disjunct_constraints(model)) == orig_num_disj_cons
+    @test length(DP._disjunctions(model)) == orig_num_disj
+    @test !haskey(object_dictionary(model), :W)
+    
+    # Store original reformulation methods
+    orig_methods = Dict(
+        :model => [DP._solution_method(model)],
+        :new_model1 => [DP._solution_method(new_model1)]
+    )
+    
+    # Solve new_model with Big-M reformulation
+    set_optimizer(new_model, HiGHS.Optimizer)
+    set_silent(new_model)
+    DP.optimize!(new_model, gdp_method = BigM())
+    @test termination_status(new_model) in (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
+    
+    # Verify reformulation methods of other models are unchanged
+    @test DP._solution_method(model) == orig_methods[:model][]
+    @test DP._solution_method(new_model1) == orig_methods[:new_model1][]
+    
+    # Verify model and new_model1 were not modified by the solve
+    @test num_variables(model) == orig_num_vars
+    @test num_variables(new_model1) == orig1_num_vars
+    @test length(DP._logical_variables(model)) == orig_num_lvars
+    @test length(DP._logical_variables(new_model1)) == orig1_num_lvars
+    @test length(DP._disjunct_constraints(model)) == orig_num_disj_cons
+    @test length(DP._disjunct_constraints(new_model1)) == orig1_num_disj_cons
+    @test length(DP._disjunctions(model)) == orig_num_disj
+    @test length(DP._disjunctions(new_model1)) == orig1_num_disj
+    
+    @test !haskey(object_dictionary(new_model1), :W)
+    
+    @test length(DP._logical_variables(new_model)) == orig_num_lvars + 2
+    @test length(DP._disjunct_constraints(new_model)) == orig_num_disj_cons + 2
+    @test length(DP._disjunctions(new_model)) == orig_num_disj + 1
+end
+
 @testset "GDP Model" begin
     test_GDPData()
     test_empty_model()
     test_non_gdp_model()
+    test_copy_model()
     test_creation_optimizer()
     test_set_optimizer()
-end
\ No newline at end of file
+    test_remapping_functions()
+end
diff --git a/test/solve.jl b/test/solve.jl
index cfd24bd..bb503bd 100644
--- a/test/solve.jl
+++ b/test/solve.jl
@@ -53,6 +53,20 @@ function test_linear_gdp_example(m, use_complements = false)
     @test value(Y[2])
     @test !value(W[1])
     @test !value(W[2])
+
+    m_copy, ref_map = JuMP.copy_model(m)
+    lv_map = DP.copy_gdp_data(m, m_copy, ref_map)
+    set_optimizer(m_copy, HiGHS.Optimizer)
+    set_optimizer_attribute(m_copy, "output_flag", false)
+    optimize!(m_copy, gdp_method = BigM())
+    @test termination_status(m_copy) == MOI.OPTIMAL
+    @test objective_value(m_copy) ≈ 11
+    @test value.(x) == value.(ref_map[x])
+    @test value.(Y[1]) == value.(lv_map[Y[1]])
+    @test value.(Y[2]) == value.(lv_map[Y[2]])
+    @test !value(W[1])
+    @test !value(W[2])
+
 end
 
 function test_quadratic_gdp_example(use_complements = false) #psplit does not work with complements