diff --git a/Project.toml b/Project.toml
index c27f5dc..ceda3d3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -12,11 +12,15 @@ Aqua = "0.8"
 JuMP = "1.18"
 Reexport = "1"
 julia = "1.6"
+Juniper = "0.9.3"
+Ipopt = "1.9.0"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
+Juniper = "2ddba703-00a4-53a7-87a5-e8b9971dde84"
 
 [targets]
-test = ["Aqua", "HiGHS", "Test"]
+test = ["Aqua", "HiGHS", "Test", "Juniper", "Ipopt"]
diff --git a/README.md b/README.md
index 09ffc22..5d68bae 100644
--- a/README.md
+++ b/README.md
@@ -177,6 +177,14 @@ The following reformulation methods are currently supported:
     - `optimizer`: Optimizer to use when solving subproblems to determine M values. This is a required value.
     - `default_M`: Default big-M value to use if no big-M is specified for a logical variable (1e9).
 
+5. [P-Split](https://arxiv.org/abs/2202.05198): This method reformulates each disjunct constraint into P constraints, each with a partitioned group defined by the user. This method requires that terms in the constraint be convex additively seperable with respect to each variable. The `PSplit` struct is created with the following required arguments:
+
+    - `partition`: Partition of the variables to be split. All variables must be in exactly one partition. (e.g., The variables `x[1:4]` can be partitioned into two groups ` partition = [[x[1], x[2]], [x[3], x[4]]]`)
+    - `PSplit(n_parts, model)`: Automatically partition all variables in the model into `n_parts` groups
+
+    All variables must be included in exactly one partition. For manual partitioning, ensure each variable appears in exactly one group. For automatic partitioning, variables are divided as evenly as possible among the specified number of partitions.
+
+
 ## Release Notes
 
 Prior to `v0.4.0`, the package did not leverage the JuMP extension capabilities and was not as robust. For these earlier releases, refer to [Perez, Joshi, and Grossmann, 2023](https://arxiv.org/abs/2304.10492v1) and the following [JuliaCon 2022 Talk](https://www.youtube.com/watch?v=AMIrgTTfUkI).
@@ -192,7 +200,7 @@ using HiGHS
 m = GDPModel(HiGHS.Optimizer)
 @variable(m, 0 ≤ x[1:2] ≤ 20)
 @variable(m, Y[1:2], Logical)
-@constraint(m, [i = 1:2], [2,5][i] ≤ x[i] ≤ [6,9][i], Disjunct(Y[1]))
+@constraint(m, [i = 1:2], [2,5][i] ≤ x[i] ≤ [6,9][i], Disjunct(Y[1])) 
 @constraint(m, [i = 1:2], [8,10][i] ≤ x[i] ≤ [11,15][i], Disjunct(Y[2]))
 @disjunction(m, Y)
 @objective(m, Max, sum(x))
diff --git a/src/DisjunctiveProgramming.jl b/src/DisjunctiveProgramming.jl
index 5e3b162..3bf91c8 100644
--- a/src/DisjunctiveProgramming.jl
+++ b/src/DisjunctiveProgramming.jl
@@ -25,6 +25,7 @@ include("mbm.jl")
 include("indicator.jl")
 include("print.jl")
 include("utilities.jl")
+include("psplit.jl")
 
 # Define additional stuff that should not be exported
 const _EXCLUDE_SYMBOLS = [Symbol(@__MODULE__), :eval, :include]
diff --git a/src/datatypes.jl b/src/datatypes.jl
index 2248018..d46f4e0 100644
--- a/src/datatypes.jl
+++ b/src/datatypes.jl
@@ -431,6 +431,72 @@ mutable struct _Hull{V <: JuMP.AbstractVariableRef, T} <: AbstractReformulationM
     end
 end
 
+"""
+    PSplit <: AbstractReformulationMethod
+
+A type for using the P-split reformulation approach for disjunctive constraints.
+This method partitions variables into groups and handles each group separately.
+
+# Constructors
+- `PSplit(partition::Vector{Vector{V}})`: Create a PSplit with the given 
+partition of variables
+- `PSplit(n_parts::Int, model::JuMP.AbstractModel)`: Automatically partition 
+model variables into `n_parts` groups
+
+# Fields
+- `partition::Vector{Vector{V}}`: The partition of variables, where each inner 
+vector represents a group of variables that will be handled together
+"""
+struct PSplit{V <: JuMP.AbstractVariableRef} <: AbstractReformulationMethod
+    partition::Vector{Vector{V}}
+
+    function PSplit(partition::Vector{Vector{V}}) where 
+        {V <: JuMP.AbstractVariableRef}
+        new{V}(partition)
+    end
+
+    function PSplit(n_parts::Int, model::JuMP.AbstractModel)
+        n_parts > 0 || error("Number of partitions must be 
+        positive, got $n_parts")
+        variables = collect(JuMP.all_variables(model))
+        n_vars = length(variables)
+        
+        n_parts = min(n_parts, n_vars)
+        n_parts > 0 || error("No variables found in the model")
+        
+        base_size = n_vars ÷ n_parts
+        remaining = n_vars % n_parts
+        
+        partition = Vector{Vector{eltype(variables)}}()
+        start_idx = 1
+        
+        for i in 1:n_parts
+            part_size = i <= remaining ? base_size + 1 : base_size
+            end_idx = start_idx + part_size - 1
+            push!(partition, variables[start_idx:end_idx])
+            start_idx = end_idx + 1
+        end
+        
+        return PSplit(partition)
+    end
+end
+
+# temp struct to store variable disaggregations (reset for each disjunction)
+mutable struct _PSplit{V <: JuMP.AbstractVariableRef, M <: JuMP.AbstractModel, T} <: AbstractReformulationMethod
+    partition::Vector{Vector{V}}
+    sum_constraints::Dict{LogicalVariableRef{M}, Vector{<:AbstractConstraint}}
+    hull::_Hull{V, T}
+    function _PSplit(method::PSplit{V}, model::M) where 
+        {V <: JuMP.AbstractVariableRef, M <: JuMP.AbstractModel}
+        T = JuMP.value_type(M)
+        new{V, M, T}(
+            method.partition, 
+            Dict{LogicalVariableRef{M}, Vector{<:AbstractConstraint}}(), 
+            _Hull(Hull(), Set{V}())
+        )
+    end
+end
+
 """
     Indicator <: AbstractReformulationMethod
 
@@ -539,4 +605,4 @@ function VariableProperties(vref::JuMP.GenericVariableRef{T}) where T
     name = JuMP.name(vref)
     set = JuMP.is_variable_in_set(vref) ? JuMP.moi_set(JuMP.constraint_object(JuMP.VariableInSetRef(vref))) : nothing
     return VariableProperties(info, name, set, nothing)
-end
\ No newline at end of file
+end
diff --git a/src/hull.jl b/src/hull.jl
index cba5efd..8ad5209 100644
--- a/src/hull.jl
+++ b/src/hull.jl
@@ -22,14 +22,14 @@ function _disaggregate_variables(
     end
 end
 function _disaggregate_variable(
-    model::JuMP.AbstractModel, 
+    model::M, 
     lvref::LogicalVariableRef, 
     vref::JuMP.AbstractVariableRef, 
     method::_Hull
-    )
+    ) where {M <: JuMP.AbstractModel}
     #create disaggregated vref
     lb, ub = variable_bound_info(vref)
-    T = JuMP.value_type(typeof(model))
+    T = JuMP.value_type(M)
     info = JuMP.VariableInfo(
         true,      # has_lb = true
         lb,        # lower_bound = lb
diff --git a/src/psplit.jl b/src/psplit.jl
new file mode 100644
index 0000000..dea6789
--- /dev/null
+++ b/src/psplit.jl
@@ -0,0 +1,492 @@
+################################################################################
+#                              BUILD PARTITIONED EXPRESSION
+################################################################################
+
+function _build_partitioned_expression(
+    expr::T,
+    partition_variables::Vector{<:JuMP.AbstractVariableRef}
+) where {T <: JuMP.GenericAffExpr}
+    constant = JuMP.constant(expr)
+    new_affexpr = zero(T)
+    for var in partition_variables
+        JuMP.add_to_expression!(new_affexpr, JuMP.coefficient(expr, var), var) 
+    end
+    return new_affexpr, constant
+end
+
+function _build_partitioned_expression(
+    expr::T,
+    partition_variables::Vector{<:JuMP.AbstractVariableRef}
+) where {T <: JuMP.GenericQuadExpr}
+    new_quadexpr = zero(T)
+    constant = JuMP.constant(expr)
+    for var in partition_variables
+        for (pair, coeff) in expr.terms
+            if pair.a == var && pair.b == var
+                JuMP.add_to_expression!(new_quadexpr, coeff, var, var)
+            elseif pair.a == var || pair.b == var
+                error("Quadratic expression contains 
+                bilinear term ($(pair.a), $(pair.b))")
+            end
+            
+        end
+    end
+    new_aff, _ = _build_partitioned_expression(expr.aff, partition_variables)
+    JuMP.add_to_expression!(new_quadexpr, new_aff)
+    return new_quadexpr, constant
+end
+
+function _build_partitioned_expression(
+    expr::T,
+    partition_variables::Vector{<:JuMP.AbstractVariableRef}
+) where {T <: JuMP.AbstractVariableRef}
+
+    if expr in partition_variables
+        return expr, zero(T)
+    else
+        return zero(T), zero(T)
+    end
+end
+
+function _build_partitioned_expression(
+    expr::T,
+    partition_variables::Vector{<:JuMP.AbstractVariableRef}
+) where {T <: Number}
+    return expr, zero(T)
+end
+
+function _build_partitioned_expression(
+    expr::T,
+    partition_variables::Vector{<:JuMP.AbstractVariableRef}
+) where {T <: JuMP.GenericNonlinearExpr}
+    error("P-Split does not currently support nonlinear expressions $(expr)")
+end
+
+
+################################################################################
+#                              BOUND AUXILIARY
+################################################################################
+function _bound_auxiliary(
+    model::M,
+    v::JuMP.AbstractVariableRef,
+    func::JuMP.GenericAffExpr{T,V},
+    method::PSplit
+) where {M <: JuMP.AbstractModel, T, V}
+
+    lower_bound = has_lower_bound(v) ? lower_bound(v) : zero(T)
+    upper_bound = has_upper_bound(v) ? upper_bound(v) : zero(T)
+
+    for (var, coeff) in func.terms
+        if var != v
+            JuMP.is_binary(var) && continue
+            var_lb, var_ub = variable_bound_info(var)
+            if coeff > 0.0
+                lower_bound += coeff * var_lb
+                upper_bound += coeff * var_ub
+            else
+                lower_bound += coeff * var_ub
+                upper_bound += coeff * var_lb
+            end
+        end
+    end
+    JuMP.set_lower_bound(v, lower_bound)
+    JuMP.set_upper_bound(v, upper_bound)
+    _variable_bounds(model)[v] = set_variable_bound_info(v, method)
+end    
+
+function _bound_auxiliary(
+    model::M,
+    v::JuMP.AbstractVariableRef,
+    func::Number,
+    method::PSplit
+) where {M <: JuMP.AbstractModel}
+    JuMP.set_lower_bound(v, func)
+    JuMP.set_upper_bound(v, func)
+    _variable_bounds(model)[v] = set_variable_bound_info(v, method)
+    return
+end
+
+function _bound_auxiliary(
+    model::M,
+    v::JuMP.AbstractVariableRef,
+    func::JuMP.GenericQuadExpr{T,V},
+    method::PSplit
+) where {M <: JuMP.AbstractModel, T, V}
+    
+    # Handle linear terms
+    _bound_auxiliary(model, v, func.aff, method)
+    lower_bound = JuMP.lower_bound(v)
+    upper_bound = JuMP.upper_bound(v)
+    
+    # Handle quadratic terms
+    for (vars, coeff) in func.terms
+        var = vars.a 
+        if var != v
+            JuMP.is_binary(var) && continue
+            lb, ub = variable_bound_info(var)
+            
+            # For x^2 terms
+            sq_min = min(lb^2, ub^2, zero(T))
+            sq_max = max(lb^2, ub^2, zero(T))
+            
+            if coeff > 0.0
+                lower_bound += coeff * sq_min
+                upper_bound += coeff * sq_max
+            else
+                lower_bound += coeff * sq_max
+                upper_bound += coeff * sq_min
+            end
+        end
+    end
+    
+    # Add constant term
+    const_term = func.aff.constant
+    # lower_bound += const_term
+    # upper_bound += const_term
+    
+    JuMP.set_lower_bound(v, lower_bound)
+    JuMP.set_upper_bound(v, upper_bound)
+    _variable_bounds(model)[v] = set_variable_bound_info(v, method)
+end
+
+function _bound_auxiliary(
+    model::M,
+    v::JuMP.AbstractVariableRef,
+    func::JuMP.AbstractVariableRef,
+    method::PSplit
+) where {M <: JuMP.AbstractModel}
+    T = JuMP.value_type(M)
+    lower_bound = zero(T)
+    upper_bound = zero(T)
+    if func != v
+        lower_bound = variable_bound_info(func)[1]
+        upper_bound = variable_bound_info(func)[2]
+        JuMP.set_lower_bound(v, lower_bound)
+        JuMP.set_upper_bound(v, upper_bound)
+    else
+        JuMP.set_lower_bound(v,lower_bound)
+        JuMP.set_upper_bound(v,upper_bound)
+    end
+    _variable_bounds(model)[v] = set_variable_bound_info(v, method)
+end
+
+requires_variable_bound_info(method::Union{PSplit, _PSplit}) = true
+
+function set_variable_bound_info(
+    vref::JuMP.AbstractVariableRef, 
+    ::Union{PSplit, _PSplit}
+    )
+    if !has_lower_bound(vref) || !has_upper_bound(vref)
+        error("Variable $vref must have both lower and upper bounds defined when
+         using the PSplit reformulation."
+         )
+    else
+        lb = min(0, lower_bound(vref))
+        ub = max(0, upper_bound(vref))
+    end
+    return lb, ub
+end
+
+################################################################################
+#                              REFORMULATE DISJUNCT
+################################################################################
+
+function reformulate_disjunction(
+    model::JuMP.AbstractModel, 
+    disj::Disjunction, 
+    method::PSplit{V}
+) where {V <: JuMP.AbstractVariableRef}
+    ref_cons = Vector{JuMP.AbstractConstraint}() #store reformulated constraints
+    disj_vrefs = _get_disjunction_variables(model, disj)
+    sum_constraints = Dict{LogicalVariableRef, Vector{<:JuMP.AbstractConstraint}}()
+    aux_vars = Set{V}()
+    for d in disj.indicators
+        partitioned_constraints, sum_constraints[d], vars = _partition_disjunct(model, d, method)
+        append!(ref_cons, partitioned_constraints)
+        union!(aux_vars, vars)
+    end
+    psplit = _PSplit(method, model)
+    psplit.hull = _Hull(Hull(), union(disj_vrefs, aux_vars))
+    psplit.sum_constraints = sum_constraints
+    for d in disj.indicators
+        bvref = binary_variable(d)
+        for vref in disj_vrefs
+            push!(psplit.hull.disjunction_variables[vref], vref)
+            psplit.hull.disjunct_variables[vref, bvref] = vref
+        end
+        _disaggregate_variables(model, d, aux_vars, psplit.hull)
+        _reformulate_disjunct(model, ref_cons, d, psplit)
+    end
+    for vref in aux_vars
+        _aggregate_variable(model, ref_cons, vref, psplit.hull)
+    end
+    return ref_cons
+end
+
+function reformulate_disjunction(
+    model::JuMP.AbstractModel, 
+    disj::Disjunction, 
+    method::_PSplit
+)
+    return reformulate_disjunction(model, disj, PSplit(method.partition))
+end
+
+function _reformulate_disjunct(
+    model::JuMP.AbstractModel, 
+    ref_cons::Vector{JuMP.AbstractConstraint}, 
+    lvref::LogicalVariableRef, 
+    method::_PSplit
+    )
+    #reformulate each constraint and add to the model
+    bvref = binary_variable(lvref)
+    haskey(method.sum_constraints, lvref) || return
+    constraints = method.sum_constraints[lvref]
+    for con in constraints
+        append!(ref_cons, reformulate_disjunct_constraint(model, con, bvref, method.hull))
+    end
+    return
+end
+
+function _partition_disjunct(
+    model::M, 
+    lvref::LogicalVariableRef, 
+    method::PSplit
+) where {M <: JuMP.AbstractModel}
+    !haskey(_indicator_to_constraints(model), lvref) && return #skip if disjunct is empty
+    
+    partitioned_constraints = Vector{AbstractConstraint}()
+    sum_constraints = Vector{AbstractConstraint}()
+    aux_vars = Set{JuMP.AbstractVariableRef}()
+    for cref in _indicator_to_constraints(model)[lvref] 
+        con = JuMP.constraint_object(cref)
+        if !(con isa Disjunction)
+            part_con, sum_con, new_aux_vars = _build_partitioned_constraint(model, con, method)
+            append!(partitioned_constraints, part_con)
+            append!(sum_constraints, sum_con)
+            union!(aux_vars, new_aux_vars)   
+        end
+    end
+    return partitioned_constraints, sum_constraints, aux_vars
+end
+
+#################################################################################
+#                              BUILD PARTITIONED CONSTRAINT
+#################################################################################
+function _build_partitioned_constraint(
+    model::M,
+    con::JuMP.ScalarConstraint{T, S},
+    method::PSplit
+) where {M <: JuMP.AbstractModel, T, S <: _MOI.LessThan}
+    val_type = JuMP.value_type(M)
+    p = length(method.partition)
+    v = [@variable(model, base_name = "v_$(hash(con))_$(i)") for i in 1:p]
+    _, constant = _build_partitioned_expression(con.func, method.partition[p])
+    part_con = Vector{JuMP.AbstractConstraint}(undef, p)
+    for i in 1:p
+        func, _ = _build_partitioned_expression(con.func, method.partition[i])
+        part_con[i] = JuMP.build_constraint(error, func - v[i], 
+        MOI.LessThan(zero(val_type))
+        )
+        _bound_auxiliary(model, v[i], func, method)
+    end
+    sum_con = JuMP.build_constraint(error, sum(v[i] for i in 1:p) 
+    ,MOI.LessThan(con.set.upper - constant)
+    )
+
+    return part_con, [sum_con], v
+end
+
+function _build_partitioned_constraint(
+    model::M,
+    con::JuMP.ScalarConstraint{T, S},
+    method::PSplit
+) where {M <: JuMP.AbstractModel, T, S <: _MOI.GreaterThan}
+    val_type = JuMP.value_type(M)
+    p = length(method.partition)
+    part_con = Vector{JuMP.AbstractConstraint}(undef, p)
+    v = [@variable(model, base_name = "v_$(hash(con))_$(i)") for i in 1:p]
+    _, constant = _build_partitioned_expression(con.func, method.partition[p])
+    for i in 1:p
+        func, _ = _build_partitioned_expression(con.func, method.partition[i])
+        part_con[i] = JuMP.build_constraint(error, -func - v[i], 
+        MOI.LessThan(zero(val_type))
+        )
+        _bound_auxiliary(model, v[i], -func, method)
+    end
+    sum_con = JuMP.build_constraint(error, sum(v[i] for i in 1:p) 
+    , MOI.LessThan(-con.set.lower + constant)
+    )
+    return part_con, [sum_con], v
+end
+
+function _build_partitioned_constraint(
+    model::M,
+    con::JuMP.ScalarConstraint{T, S},
+    method::PSplit
+) where {M <: JuMP.AbstractModel, T, S <: Union{_MOI.Interval, _MOI.EqualTo}}
+    val_type = JuMP.value_type(M)
+    p = length(method.partition)
+    part_con_lt = Vector{JuMP.AbstractConstraint}(undef, p)
+    part_con_gt = Vector{JuMP.AbstractConstraint}(undef, p)
+    #let [_, 1] be the upper bound and [_, 2] be the lower bound
+    _, constant = _build_partitioned_expression(con.func, method.partition[p]) 
+    v = [@variable(
+        model, 
+        base_name = "v_$(hash(con))_$(i)_$(j)"
+        ) for i in 1:p, j in 1:2]
+    for i in 1:p
+        func, _= _build_partitioned_expression(con.func, method.partition[i])
+        part_con_lt[i] = JuMP.build_constraint(error, 
+        func - v[i,1], MOI.LessThan(zero(val_type))
+        )
+        part_con_gt[i] = JuMP.build_constraint(error, 
+        -func - v[i,2], MOI.LessThan(zero(val_type))
+        )
+        _bound_auxiliary(model, v[i,1], func, method)
+        _bound_auxiliary(model, v[i,2], -func, method)
+    end
+    set_values = _set_values(con.set)
+    sum_con_lt = JuMP.build_constraint(error, sum(v[i,1] for i in 1:p), 
+    MOI.LessThan((set_values[2] - constant))
+    )
+    sum_con_gt = JuMP.build_constraint(error, sum(v[i,2] for i in 1:p), 
+    MOI.LessThan(-set_values[1] + constant)
+     )
+    return vcat(part_con_lt, part_con_gt), [sum_con_lt, sum_con_gt], vec(v)
+end
+function _build_partitioned_constraint(
+    model::M,
+    con::JuMP.VectorConstraint{T, S, R},
+    method::PSplit
+) where {M <: JuMP.AbstractModel, T, S <: _MOI.Nonpositives, R}
+    p = length(method.partition)
+    d = con.set.dimension
+    v = [@variable(
+        model, 
+        base_name = "v_$(hash(con))_$(i)_$(j)"
+        ) for i in 1:p, j in 1:d]
+    part_con = Vector{JuMP.AbstractConstraint}(undef, p)
+    constants = Vector{Number}(undef, d)
+    for i in 1:p
+        part_expr = [_build_partitioned_expression(con.func[j],
+        method.partition[i]) for j in 1:d
+        ]
+        func = JuMP.@expression(model, [j = 1:d], part_expr[j][1])
+        constants .= [part_expr[j][2] for j in 1:d]
+        part_con[i] = JuMP.build_constraint(error, 
+        func - v[i,:], _MOI.Nonpositives(d)
+        )
+        for j in 1:d
+            _bound_auxiliary(model, v[i,j], func[j], method)
+        end
+    end
+    new_func = JuMP.@expression(model,[j = 1:d], 
+    sum(v[i,j] for i in 1:p) + constants[j]
+    )
+    sum_con = JuMP.build_constraint(error, new_func, _MOI.Nonpositives(d))
+    return part_con, [sum_con], vec(v)
+end
+
+function _build_partitioned_constraint(
+    model::M,
+    con::JuMP.VectorConstraint{T, S, R},
+    method::PSplit
+) where {M <: JuMP.AbstractModel, T, S <: _MOI.Nonnegatives, R}
+    p = length(method.partition)
+    d = con.set.dimension
+    v = [@variable(
+        model, 
+        base_name = "v_$(hash(con))_$(i)_$(j)"
+        ) for i in 1:p, j in 1:d]
+    part_con = Vector{JuMP.AbstractConstraint}(undef, p)
+    constants = Vector{Number}(undef, d)
+    for i in 1:p
+        part_expr = [
+            _build_partitioned_expression(con.func[j], method.partition[i]) 
+            for j in 1:d
+        ]
+        func = JuMP.@expression(model, [j = 1:d], -part_expr[j][1])
+        constants .= [-part_expr[j][2] for j in 1:d]
+        part_con[i] = JuMP.build_constraint(error, func - v[i,:], _MOI.Nonpositives(d))
+        for j in 1:d
+            _bound_auxiliary(model, v[i,j], func[j], method)
+        end
+    end
+    new_func = JuMP.@expression(model,[j = 1:d], 
+    sum(v[i,j] for i in 1:p) + constants[j]
+    )
+    sum_con = JuMP.build_constraint(error,new_func,_MOI.Nonpositives(d))
+    return part_con, [sum_con], vec(v)
+end
+
+function _build_partitioned_constraint(
+    model::M,
+    con::JuMP.VectorConstraint{T, S, R},
+    method::PSplit
+) where {M <: JuMP.AbstractModel, T, S <: _MOI.Zeros, R}
+    p = length(method.partition)
+    d = con.set.dimension
+    part_con_np = Vector{JuMP.AbstractConstraint}(undef, p)  # nonpositive (≤ 0)
+    part_con_nn = Vector{JuMP.AbstractConstraint}(undef, p)  # nonnegative (≥ 0)
+    v = [@variable(
+        model, 
+        base_name = "v_$(hash(con))_$(i)_$(j)_$(k)"
+        ) for i in 1:p, j in 1:d, k in 1:2
+    ]
+    constants = Vector{Number}(undef, d)
+    for i in 1:p
+        part_expr = [
+            _build_partitioned_expression(con.func[j], method.partition[i]) 
+            for j in 1:d
+        ]
+        func = JuMP.@expression(model, [j = 1:d], part_expr[j][1])        
+        constants .= [part_expr[j][2] for j in 1:d]
+        part_con_np[i] = JuMP.build_constraint(error, 
+        func - v[i,:,1], _MOI.Nonpositives(d)
+        )
+        part_con_nn[i] = JuMP.build_constraint(error, 
+        -func - v[i,:,2], _MOI.Nonpositives(d)
+        )
+        for j in 1:d
+            _bound_auxiliary(model, v[i,j,1], func[j], method)
+            _bound_auxiliary(model, v[i,j,2], -func[j], method)
+        end
+    end
+    new_func_np = JuMP.@expression(model,[j = 1:d], 
+    sum(v[i,j,1] for i in 1:p) + constants[j]
+    )
+    new_func_nn = JuMP.@expression(model,[j = 1:d], 
+    -sum(v[i,j,2] for i in 1:p) - constants[j]
+    )
+    sum_con_np = JuMP.build_constraint(error, 
+    new_func_np, _MOI.Nonpositives(d)
+    )
+    sum_con_nn = JuMP.build_constraint(error, 
+    new_func_nn, _MOI.Nonpositives(d)
+    )
+    return vcat(part_con_np, part_con_nn), vcat(sum_con_np, sum_con_nn), vec(v)
+end
+
+################################################################################
+#                          FALLBACK WARNING DISPATCHES
+################################################################################
+
+# Generic fallback for _build_partitioned_expression
+function _build_partitioned_expression(
+    expr::F,
+    ::Vector{<:JuMP.AbstractVariableRef}
+) where F
+    error("PSplit: _build_partitioned_expression not implemented 
+    for expression type $F.")
+end
+
+# Generic fallback for _bound_auxiliary
+function _bound_auxiliary(
+    ::JuMP.AbstractModel,
+    v::JuMP.AbstractVariableRef,
+    func::F,
+    ::PSplit
+) where F
+    error("PSplit: _bound_auxiliary not implemented for function 
+    type $F.")
+end
\ No newline at end of file
diff --git a/test/constraints/psplit.jl b/test/constraints/psplit.jl
new file mode 100644
index 0000000..cc07cb5
--- /dev/null
+++ b/test/constraints/psplit.jl
@@ -0,0 +1,292 @@
+function test_psplit()
+    model = GDPModel()
+    @variable(model, x[1:4])
+    method = PSplit([[x[1], x[2]], [x[3], x[4]]])
+    @test method.partition == [[x[1], x[2]], [x[3], x[4]]]
+
+    method = PSplit(3, model)
+    @test length(method.partition) == 3
+    @test method.partition[1] == [x[1], x[2]]
+    @test method.partition[2] == [x[3]]
+    @test method.partition[3] == [x[4]]
+
+end
+
+function test_build_partitioned_expression()
+    model = JuMP.Model()
+    @variable(model, x[1:4])
+    partition_variables = [x[1], x[4]]
+    
+    test_cases = [
+        (expr = 1.0 * x[1] - 2.0 * x[2], expected = (1.0 * x[1], 0.0)),
+        (expr = x[1] * x[1] + 2.0 * x[1] * x[1] + 3.0 * x[2] * x[2], 
+        expected = (3 * x[1] * x[1], 0.0)),
+        (expr = x[3], expected = (0, 0)),
+        (expr = x[4], expected = (x[4], 0.0)),
+        (expr = 4.0, expected = (4.0, 0.0))
+    ]
+
+    for (i, test_case) in enumerate(test_cases)
+        result = DP._build_partitioned_expression(test_case.expr, partition_variables)
+        @test result == test_case.expected
+    end
+
+    @test_throws ErrorException DP._build_partitioned_expression(
+        x[1] * x[2], 
+        partition_variables
+    )
+    @test_throws ErrorException DP._build_partitioned_expression(
+        exp(x[1]), 
+        partition_variables
+    )
+    @test_throws ErrorException DP._build_partitioned_expression(
+        "Bad Input", 
+        partition_variables
+    )
+end
+
+
+
+
+
+function test_bound_auxiliary()
+    model = GDPModel()
+    @variable(model, 0 <= x[1:4] <= 3)
+    @variable(model, v[1:6])
+    method = PSplit([[x[1], x[2]], [x[3], x[4]]])
+    affexpr = 1.0 * x[1] - 2.0 * x[2]
+    quadexpr = x[1] * x[1] + 2.0 * x[1] * x[1] - 3.0 * x[2] * x[2] + x[1] - x[2]
+    var = x[3]
+    num = 4.0
+    nl = exp(x[1])
+    for i in 1:4
+        DP._variable_bounds(model)[x[i]] = DP.set_variable_bound_info(x[i], method)
+    end
+
+    DP._bound_auxiliary(model, v[2], quadexpr, method)
+    DP._bound_auxiliary(model, v[3], num, method)
+    DP._bound_auxiliary(model, v[4], var, method)
+    DP._bound_auxiliary(model, v[5], affexpr, method)
+    DP._bound_auxiliary(model, v[6], v[6], method)
+
+    @test JuMP.lower_bound(v[6]) == 0
+    @test JuMP.upper_bound(v[6]) == 0
+    @test JuMP.lower_bound(v[5]) == -6
+    @test JuMP.upper_bound(v[5]) == 3
+    @test JuMP.lower_bound(v[4]) == 0
+    @test JuMP.upper_bound(v[4]) == 3
+    @test JuMP.lower_bound(v[2]) == -30
+    @test JuMP.upper_bound(v[2]) == 30
+
+    @test JuMP.upper_bound(v[3]) == 4
+    @test JuMP.lower_bound(v[3]) == 4
+    @test_throws ErrorException DP._bound_auxiliary(model, v[3], nl, method)
+end
+
+
+function test_build_partitioned_constraint()
+    model = GDPModel()
+    @variable(model, 0 <= x[1:4] <= 3)
+    @variable(model, y[1:2], Bin)
+    method = PSplit([[x[1], x[2]], [x[3], x[4]]])
+    for i in 1:4
+        DP._variable_bounds(model)[x[i]] = DP.set_variable_bound_info(
+            x[i], 
+            method
+        )
+    end
+
+    lt = JuMP.constraint_object(@constraint(model, x[1] + x[3] <= 1))
+    gt = JuMP.constraint_object(@constraint(model, x[2] + x[4] >= 1))
+    eq = JuMP.constraint_object(@constraint(model, x[3] + x[4] == 1))
+    interval = JuMP.constraint_object(
+        @constraint(model, 0 <= x[1] + x[2] + x[3] <= 0.5)
+    )
+    nn = JuMP.constraint_object(@constraint(model, x .- 1 >= 0))
+    np = JuMP.constraint_object(@constraint(model, -x .+ 1 <= 0))
+    zeros = JuMP.constraint_object(@constraint(model, 5x .- 1 == 0))
+    
+    ref_lt, lt_sum, lt_vars = DP._build_partitioned_constraint(model, lt, method)
+    ref_gt, gt_sum, gt_vars = DP._build_partitioned_constraint(model, gt, method)
+    ref_eq, eq_sum, eq_vars = DP._build_partitioned_constraint(model, eq, method)
+    ref_interval, interval_sum, interval_vars = DP._build_partitioned_constraint(
+        model, interval, method
+    )
+    ref_nn, nn_sum, nn_vars = DP._build_partitioned_constraint(model, nn, method)
+    ref_np, np_sum, np_vars = DP._build_partitioned_constraint(model, np, method)
+    ref_zeros, zeros_sum, zeros_vars = DP._build_partitioned_constraint(model, zeros, method)
+
+    @test ref_lt[1].func == x[1] - 
+        variable_by_name(model, "v_$(hash(lt))_1")
+    @test ref_lt[2].func == x[3] - 
+        variable_by_name(model, "v_$(hash(lt))_2")
+    @test lt_sum[1].func == 
+        variable_by_name(model, "v_$(hash(lt))_1") + 
+        variable_by_name(model, "v_$(hash(lt))_2")
+    @test ref_gt[1].func == -x[2] - 
+        variable_by_name(model, "v_$(hash(gt))_1")
+    @test ref_gt[2].func == -x[4] - 
+        variable_by_name(model, "v_$(hash(gt))_2")
+    @test gt_sum[1].func == 
+        variable_by_name(model, "v_$(hash(gt))_1") + 
+        variable_by_name(model, "v_$(hash(gt))_2")
+    @test ref_eq[1].func == 
+        -variable_by_name(model, "v_$(hash(eq))_1_1")
+    @test ref_eq[2].func == x[3] + x[4] - 
+        variable_by_name(model, "v_$(hash(eq))_2_1")
+    @test ref_eq[4].func == -x[3] - x[4] - 
+        variable_by_name(model, "v_$(hash(eq))_2_2")
+    @test eq_sum[1].func == 
+        variable_by_name(model, "v_$(hash(eq))_1_1") + 
+        variable_by_name(model, "v_$(hash(eq))_2_1")
+    @test eq_sum[2].func == 
+        variable_by_name(model, "v_$(hash(eq))_1_2") + 
+        variable_by_name(model, "v_$(hash(eq))_2_2")
+    @test ref_interval[1].func == x[1] + x[2] - 
+        variable_by_name(model, "v_$(hash(interval))_1_1")
+    @test ref_interval[2].func == x[3] - 
+        variable_by_name(model, "v_$(hash(interval))_2_1")
+    @test ref_interval[4].func == -x[3] - 
+        variable_by_name(model, "v_$(hash(interval))_2_2")
+    @test interval_sum[1].func == 
+        variable_by_name(model, "v_$(hash(interval))_1_1") + 
+        variable_by_name(model, "v_$(hash(interval))_2_1")
+
+    @test ref_nn[1].func == [
+        -x[1] - variable_by_name(model, "v_$(hash(nn))_1_1"),
+        -x[2] - variable_by_name(model, "v_$(hash(nn))_1_2"),
+        -variable_by_name(model, "v_$(hash(nn))_1_3"),
+        -variable_by_name(model, "v_$(hash(nn))_1_4")
+    ]
+
+    @test ref_nn[2].func == [
+        -variable_by_name(model, "v_$(hash(nn))_2_1"),
+        -variable_by_name(model, "v_$(hash(nn))_2_2"),
+        -x[3] - variable_by_name(model, "v_$(hash(nn))_2_3"),
+        -x[4] - variable_by_name(model, "v_$(hash(nn))_2_4")
+    ]
+
+    @test nn_sum[1].func == [
+    variable_by_name(model, "v_$(hash(nn))_1_1") + 
+    variable_by_name(model, "v_$(hash(nn))_2_1") + 1,
+    variable_by_name(model, "v_$(hash(nn))_1_2") + 
+    variable_by_name(model, "v_$(hash(nn))_2_2") + 1,
+    variable_by_name(model, "v_$(hash(nn))_1_3") + 
+    variable_by_name(model, "v_$(hash(nn))_2_3") + 1,
+    variable_by_name(model, "v_$(hash(nn))_1_4") + 
+    variable_by_name(model, "v_$(hash(nn))_2_4") + 1
+]
+
+
+    @test ref_np[1].func == [
+        -x[1] - variable_by_name(model, "v_$(hash(np))_1_1"),
+        -x[2] - variable_by_name(model, "v_$(hash(np))_1_2"),
+        -variable_by_name(model, "v_$(hash(np))_1_3"),
+        -variable_by_name(model, "v_$(hash(np))_1_4")
+    ]
+
+    @test ref_np[2].func == [
+        -variable_by_name(model, "v_$(hash(np))_2_1"),
+        -variable_by_name(model, "v_$(hash(np))_2_2"),
+        -x[3] - variable_by_name(model, "v_$(hash(np))_2_3"),
+        -x[4] - variable_by_name(model, "v_$(hash(np))_2_4")
+    ]
+
+    @test np_sum[1].func == [
+    variable_by_name(model, "v_$(hash(np))_1_1") + 
+    variable_by_name(model, "v_$(hash(np))_2_1") + 1,
+    variable_by_name(model, "v_$(hash(np))_1_2") + 
+    variable_by_name(model, "v_$(hash(np))_2_2") + 1,
+    variable_by_name(model, "v_$(hash(np))_1_3") + 
+    variable_by_name(model, "v_$(hash(np))_2_3") + 1,
+    variable_by_name(model, "v_$(hash(np))_1_4") + 
+    variable_by_name(model, "v_$(hash(np))_2_4") + 1
+]
+
+
+    @test ref_zeros[1].func == [
+        5x[1] - variable_by_name(model, "v_$(hash(zeros))_1_1_1"),
+        5x[2] - variable_by_name(model, "v_$(hash(zeros))_1_2_1"),
+        -variable_by_name(model, "v_$(hash(zeros))_1_3_1"),
+        -variable_by_name(model, "v_$(hash(zeros))_1_4_1")
+    ]
+
+    @test ref_zeros[2].func == [
+        -variable_by_name(model, "v_$(hash(zeros))_2_1_1"),
+        -variable_by_name(model, "v_$(hash(zeros))_2_2_1"),
+        5x[3] - variable_by_name(model, "v_$(hash(zeros))_2_3_1"),
+        5x[4] - variable_by_name(model, "v_$(hash(zeros))_2_4_1")
+    ]
+end
+
+function test_partition_disjunct()
+    model = GDPModel()
+    @variable(model, 0<= x[1:4] <= 10)
+    @variable(model, Y[1:2], Logical)
+    method = PSplit([[x[1], x[2]], [x[3], x[4]]])
+    
+    con1 = @constraint(model, x[1] + x[2] <= 1, Disjunct(Y[1]))
+
+    for i in 1:4
+        DP._variable_bounds(model)[x[i]] = DP.set_variable_bound_info(x[i], method)
+        DP._bound_auxiliary(model, x[i], x[i], method)
+    end
+    
+    partitioned_constraints, sum_constraints, aux_vars = DP._partition_disjunct(model, Y[1], method)
+    @test partitioned_constraints[1].func == x[1] + x[2] - variable_by_name(model, "v_$(hash(JuMP.constraint_object(con1)))_1")
+    @test partitioned_constraints[2].func == -variable_by_name(model, "v_$(hash(JuMP.constraint_object(con1)))_2")
+
+    @test sum_constraints[1].func == variable_by_name(model, "v_$(hash(JuMP.constraint_object(con1)))_1") + 
+                                 variable_by_name(model, "v_$(hash(JuMP.constraint_object(con1)))_2")
+    @test sum_constraints[1].set == MOI.LessThan(1.0)
+    @test aux_vars == Set{JuMP.AbstractVariableRef}([variable_by_name(model, "v_$(hash(JuMP.constraint_object(con1)))_1"), variable_by_name(model, "v_$(hash(JuMP.constraint_object(con1)))_2")])
+end
+
+function test_reformulate_disjunction()
+    model = GDPModel()
+    @variable(model, 0 <= x[1:4] <= 1)
+    @variable(model, Y[1:2], Logical)
+    @variable(model, W[1:2], Logical)
+    @disjunction(model, [Y[1], Y[2]])
+    @disjunction(model, [W[1], W[2]])
+    method = PSplit([[x[1], x[2]], [x[3], x[4]]])
+    for i in 1:4
+        DP._variable_bounds(model)[x[i]] = DP.set_variable_bound_info(x[i], method)
+    end
+    ref_cons1 = Vector{JuMP.AbstractConstraint}()
+    ref_cons2 = Vector{JuMP.AbstractConstraint}()
+    
+    good_quad = @constraint(model, x[1]^2 + x[2]^2 <= 1, Disjunct(W[2]))
+    good_quad2 = @constraint(model, x[3]^2 + x[4]^2 <= 1, Disjunct(W[2]))
+    affexpr = @constraint(model, x[1] + x[2] <= 1, Disjunct(W[1]))
+    nonlinear_con = @constraint(model, exp(x[1]) <= 1, Disjunct(Y[1]))
+    bad_quad = @constraint(model, x[1]*x[2] <= 1, Disjunct(Y[2]))
+
+    @test_throws ErrorException DP._reformulate_disjunct(model, ref_cons1, Y[2], method)
+    @test_throws ErrorException DP._reformulate_disjunct(model, ref_cons2, Y[1], method)
+    
+    DP._reformulate_disjunct(model, ref_cons1, W[2], method)
+    DP._reformulate_disjunct(model, ref_cons2, W[1], method)
+
+    @test length(ref_cons1) == 6
+    @test length(ref_cons2) == 3
+end
+
+function test_set_variable_bound_info()
+    model = GDPModel()
+    @variable(model, x)
+    
+    @test_throws ErrorException DP.set_variable_bound_info(x, PSplit([[x]]))
+    JuMP.set_lower_bound(x, 0.0)
+    JuMP.set_upper_bound(x, 3.0)
+    @test DP.set_variable_bound_info(x, PSplit([[x]])) == (0, 3)
+end
+ 
+@testset "P-Split Reformulation" begin
+    test_psplit()
+    test_build_partitioned_expression()
+    test_set_variable_bound_info()
+    test_bound_auxiliary()
+    test_build_partitioned_constraint()
+    test_partition_disjunct()
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index c653a3f..843a95f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -16,8 +16,9 @@ include("constraints/disjunct.jl")
 include("constraints/indicator.jl")
 include("constraints/mbm.jl")
 include("constraints/bigm.jl")
+include("constraints/psplit.jl")
 include("constraints/hull.jl")
 include("constraints/fallback.jl")
 include("constraints/disjunction.jl")
 include("print.jl")
-include("solve.jl")
+include("solve.jl")
\ No newline at end of file
diff --git a/test/solve.jl b/test/solve.jl
index 9f79f8e..cfd24bd 100644
--- a/test/solve.jl
+++ b/test/solve.jl
@@ -1,5 +1,4 @@
-using HiGHS
-
+using HiGHS, Ipopt, Juniper
 function test_linear_gdp_example(m, use_complements = false)
     set_attribute(m, MOI.Silent(), true)
     @variable(m, 1 ≤ x[1:2] ≤ 9)
@@ -54,7 +53,65 @@ function test_linear_gdp_example(m, use_complements = false)
     @test value(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
+    ipopt = optimizer_with_attributes(Ipopt.Optimizer,"print_level"=>0,"sb"=>"yes")
+    optimizer = optimizer_with_attributes(Juniper.Optimizer, "nl_solver"=>ipopt)
+    m = GDPModel(optimizer)
+    set_attribute(m, MOI.Silent(), true)
+    @variable(m, 0 ≤ x[1:2] ≤ 10)
+    
+    if use_complements
+        @variable(m, Y1, Logical)
+        @variable(m, Y2, Logical, logical_complement = Y1)
+        Y = [Y1, Y2]
+    else
+        @variable(m, Y[1:2], Logical)
+    end
+    @variable(m, W[1:2], Logical)
+    
+    @objective(m, Max, sum(x))
+    
+    @constraint(m, y1_quad, x[1]^2 + x[2]^2 ≤ 16, Disjunct(Y[1]))
+    @constraint(m, w1[i=1:2], [1, 2][i] ≤ x[i] ≤ [3, 4][i], Disjunct(W[1]))
+    @constraint(m, w1_quad, x[1]^2 ≥ 2, Disjunct(W[1]))
+    
+    @constraint(m, w2[i=1:2], [2, 1][i] ≤ x[i] ≤ [4, 3][i], Disjunct(W[2]))
+    @constraint(m, w2_quad, x[1]^2 + x[2] ≤ 10, Disjunct(W[2]))
+    
+    @constraint(m, y2_quad, x[1]^2 + 2*x[2]^2 ≤ 25, Disjunct(Y[2]))
+    @constraint(m, y2[i=1:2], [3, 2][i] ≤ x[i] ≤ [5, 3][i], Disjunct(Y[2]))
+    
+    @disjunction(m, inner, [W[1], W[2]], Disjunct(Y[1]))
+    @disjunction(m, outer, [Y[1], Y[2]])
+    
+    @test optimize!(m, gdp_method = BigM()) isa Nothing
+    @test termination_status(m) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
+    @test objective_value(m) ≈ 6.1237 atol=1e-3  
+    @test value.(x) ≈ [4.0825, 2.0412] atol=1e-3 
+    @test !value(Y[1]) 
+    @test value(Y[2])
+    @test !value(W[1]) 
+    @test !value(W[2])
 
+    @test optimize!(m, gdp_method = MBM(optimizer)) isa Nothing
+    @test termination_status(m) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
+    @test objective_value(m) ≈ 6.1237 atol=1e-3  
+    @test value.(x) ≈ [4.0825, 2.0412] atol=1e-3 
+    @test !value(Y[1]) 
+    @test value(Y[2])
+    @test !value(W[1]) 
+    @test !value(W[2])
+
+    @test optimize!(m, gdp_method = PSplit(2,m)) isa Nothing
+    @test termination_status(m) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
+    @test objective_value(m) ≈ 6.1237 atol=1e-3  
+    @test value.(x) ≈ [4.0825, 2.0412] atol=1e-3 
+    @test !value(Y[1]) 
+    @test value(Y[2])
+    @test !value(W[1]) 
+    @test !value(W[2])
 end
 
 function test_generic_model(m)
@@ -84,5 +141,6 @@ end
         MOI.Utilities.UniversalFallback(MOIU.Model{Float32}()),
         eval_objective_value = false
         )
+    test_quadratic_gdp_example()
     test_generic_model(GDPModel{Float32}(mockoptimizer))
-end
\ No newline at end of file
+end