diff --git a/src/ThreePhasePowerModels.jl b/src/ThreePhasePowerModels.jl
index 5731d4152..0fb28a495 100644
--- a/src/ThreePhasePowerModels.jl
+++ b/src/ThreePhasePowerModels.jl
@@ -14,6 +14,7 @@ const LOGGER = getlogger(PowerModels)
 include("core/constraint_template.jl")
 include("core/ref.jl")
 include("core/variable.jl")
+include("core/objective.jl")
 
 include("form/acp.jl")
 include("form/dcp.jl")
@@ -31,5 +32,6 @@ include("prob/tp_opf.jl")
 include("prob/tp_opf_bf.jl")
 include("prob/tp_ots.jl")
 include("prob/tp_pf.jl")
+include("prob/tp_debug.jl")
 
 end
diff --git a/src/core/constraint_template.jl b/src/core/constraint_template.jl
index 40ef5db5f..6897f0f5b 100644
--- a/src/core/constraint_template.jl
+++ b/src/core/constraint_template.jl
@@ -1,6 +1,32 @@
 # Three-phase specific constraints
 
 
+""
+function constraint_kcl_shunt_slack(pm::GenericPowerModel, i::Int; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
+    if !haskey(con(pm, nw, cnd), :kcl_p)
+        con(pm, nw, cnd)[:kcl_p] = Dict{Int,ConstraintRef}()
+    end
+    if !haskey(con(pm, nw, cnd), :kcl_q)
+        con(pm, nw, cnd)[:kcl_q] = Dict{Int,ConstraintRef}()
+    end
+
+    bus = ref(pm, nw, :bus, i)
+    bus_arcs = ref(pm, nw, :bus_arcs, i)
+    bus_arcs_dc = ref(pm, nw, :bus_arcs_dc, i)
+    bus_gens = ref(pm, nw, :bus_gens, i)
+    bus_loads = ref(pm, nw, :bus_loads, i)
+    bus_shunts = ref(pm, nw, :bus_shunts, i)
+
+    bus_pd = Dict(k => ref(pm, nw, :load, k, "pd", cnd) for k in bus_loads)
+    bus_qd = Dict(k => ref(pm, nw, :load, k, "qd", cnd) for k in bus_loads)
+
+    bus_gs = Dict(k => ref(pm, nw, :shunt, k, "gs", cnd) for k in bus_shunts)
+    bus_bs = Dict(k => ref(pm, nw, :shunt, k, "bs", cnd) for k in bus_shunts)
+
+    constraint_kcl_shunt_slack(pm, nw, cnd, i, bus_arcs, bus_arcs_dc, bus_gens, bus_pd, bus_qd, bus_gs, bus_bs)
+end
+
+
 ""
 function constraint_tp_voltage(pm::GenericPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
     constraint_tp_voltage(pm, nw, cnd)
diff --git a/src/core/objective.jl b/src/core/objective.jl
new file mode 100644
index 000000000..9bd70e845
--- /dev/null
+++ b/src/core/objective.jl
@@ -0,0 +1,10 @@
+"a quadratic penalty for bus power slack variables"
+function objective_min_bus_power_slack(pm::GenericPowerModel)
+    return @objective(pm.model, Min, 
+        sum(
+            sum(
+                sum( var(pm, n, c, :p_slack, i)^2 + var(pm, n, c, :q_slack, i)^2 for (i,bus) in nw_ref[:bus])
+            for c in conductor_ids(pm, n))
+        for (n, nw_ref) in nws(pm))
+    )
+end
diff --git a/src/core/variable.jl b/src/core/variable.jl
index 759472c66..9cde18556 100644
--- a/src/core/variable.jl
+++ b/src/core/variable.jl
@@ -1,3 +1,28 @@
+"generates variables for both `active` and `reactive` slack at each bus"
+function variable_bus_power_slack(pm::GenericPowerModel; kwargs...)
+    variable_active_bus_power_slack(pm; kwargs...)
+    variable_reactive_bus_power_slack(pm; kwargs...)
+end
+
+""
+function variable_active_bus_power_slack(pm::GenericPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
+    var(pm, nw, cnd)[:p_slack] = @variable(pm.model,
+        [i in ids(pm, nw, :bus)], basename="$(nw)_$(cnd)_p_slack",
+        start = PMs.getval(ref(pm, nw, :bus, i), "p_slack_start", cnd)
+    )
+end
+
+""
+function variable_reactive_bus_power_slack(pm::GenericPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd)
+    var(pm, nw, cnd)[:q_slack] = @variable(pm.model,
+        [i in ids(pm, nw, :bus)], basename="$(nw)_$(cnd)_q_slack",
+        start = PMs.getval(ref(pm, nw, :bus, i), "q_slack_start", cnd)
+    )
+end
+
+
+
+
 "variable: `w[i] >= 0` for `i` in `bus`es"
 function variable_tp_voltage_magnitude_sqr(pm::GenericPowerModel; nw::Int=pm.cnw, cnd::Int=pm.ccnd, bounded=true)
     bus_cnd = [(i, c) for i in ids(pm, nw, :bus) for c in PMs.conductor_ids(pm)]
diff --git a/src/form/acp.jl b/src/form/acp.jl
index bf5c2f414..fa3678b06 100644
--- a/src/form/acp.jl
+++ b/src/form/acp.jl
@@ -5,12 +5,28 @@ function variable_tp_voltage(pm::GenericPowerModel{T}; kwargs...) where T <: PMs
     PMs.variable_voltage_magnitude(pm; kwargs...)
 end
 
-
 "do nothing, this model does not have complex voltage constraints"
 function constraint_tp_voltage(pm::GenericPowerModel{T}; kwargs...) where T <: PMs.AbstractACPForm
 end
 
 
+""
+function constraint_kcl_shunt_slack(pm::GenericPowerModel{T}, n::Int, c::Int, i::Int, bus_arcs, bus_arcs_dc, bus_gens, bus_pd, bus_qd, bus_gs, bus_bs) where T <: PMs.AbstractACPForm
+    vm = var(pm, n, c, :vm, i)
+    p_slack = var(pm, n, c, :p_slack, i)
+    q_slack = var(pm, n, c, :q_slack, i)
+    p = var(pm, n, c, :p)
+    q = var(pm, n, c, :q)
+    pg = var(pm, n, c, :pg)
+    qg = var(pm, n, c, :qg)
+    p_dc = var(pm, n, c, :p_dc)
+    q_dc = var(pm, n, c, :q_dc)
+
+    con(pm, n, c, :kcl_p)[i] = @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs))*vm^2 + p_slack)
+    con(pm, n, c, :kcl_q)[i] = @constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs))*vm^2 + q_slack)
+end
+
+
 """
 Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form)
 
diff --git a/src/form/shared.jl b/src/form/shared.jl
index 8673ec409..9136a6a96 100644
--- a/src/form/shared.jl
+++ b/src/form/shared.jl
@@ -1,6 +1,24 @@
 # Three-phase specific constraints
 
 
+""
+function constraint_kcl_shunt_slack(pm::GenericPowerModel{T}, n::Int, c::Int, i, bus_arcs, bus_arcs_dc, bus_gens, bus_pd, bus_qd, bus_gs, bus_bs) where T <: PMs.AbstractWForms
+    w    = var(pm, n, c, :w, i)
+    p_slack = var(pm, n, c, :p_slack, i)
+    q_slack = var(pm, n, c, :q_slack, i)
+    pg   = var(pm, n, c, :pg)
+    qg   = var(pm, n, c, :qg)
+    p    = var(pm, n, c, :p)
+    q    = var(pm, n, c, :q)
+    p_dc = var(pm, n, c, :p_dc)
+    q_dc = var(pm, n, c, :q_dc)
+
+    @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs))*w + p_slack)
+    @constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs))*w + q_slack)
+end
+
+
+
 """
 Creates Ohms constraints (yt post fix indicates that Y and T values are in rectangular form)
 """
diff --git a/src/prob/tp_debug.jl b/src/prob/tp_debug.jl
new file mode 100644
index 000000000..800b7f886
--- /dev/null
+++ b/src/prob/tp_debug.jl
@@ -0,0 +1,146 @@
+# These problem formulations are used to debug Three Phase datasets 
+# that do not converge using the standard formulations
+
+export run_tp_opf_pbs, run_tp_pf_pbs
+
+""
+function run_tp_opf_pbs(data::Dict{String,Any}, model_constructor, solver; kwargs...)
+    return PMs.run_generic_model(data, model_constructor, solver, post_tp_opf_pbs; multiconductor=true, solution_builder = get_pbs_solution, kwargs...)
+end
+
+
+""
+function run_tp_opf_pbs(file::String, model_constructor, solver; kwargs...)
+    data = ThreePhasePowerModels.parse_file(file)
+    return PMs.run_generic_model(data, model_constructor, solver, post_tp_opf_pbs; multiconductor=true, solution_builder = get_pbs_solution, kwargs...)
+end
+
+
+""
+function post_tp_opf_pbs(pm::GenericPowerModel)
+    for c in PMs.conductor_ids(pm)
+        variable_tp_voltage(pm, cnd=c)
+        variable_bus_power_slack(pm, cnd=c)
+        PMs.variable_generation(pm, cnd=c)
+        PMs.variable_branch_flow(pm, cnd=c)
+        PMs.variable_dcline_flow(pm, cnd=c)
+    end
+
+    for c in PMs.conductor_ids(pm)
+        constraint_tp_voltage(pm, cnd=c)
+
+        for i in ids(pm, :ref_buses)
+            constraint_tp_theta_ref(pm, i, cnd=c)
+        end
+
+        for i in ids(pm, :bus)
+            constraint_kcl_shunt_slack(pm, i, cnd=c)
+        end
+
+        for i in ids(pm, :branch)
+            constraint_ohms_tp_yt_from(pm, i, cnd=c)
+            constraint_ohms_tp_yt_to(pm, i, cnd=c)
+
+            PMs.constraint_voltage_angle_difference(pm, i, cnd=c)
+
+            PMs.constraint_thermal_limit_from(pm, i, cnd=c)
+            PMs.constraint_thermal_limit_to(pm, i, cnd=c)
+        end
+
+        for i in ids(pm, :dcline)
+            PMs.constraint_dcline(pm, i, cnd=c)
+        end
+    end
+
+    objective_min_bus_power_slack(pm)
+end
+
+
+
+""
+function run_tp_pf_pbs(data::Dict{String,Any}, model_constructor, solver; kwargs...)
+    return PMs.run_generic_model(data, model_constructor, solver, post_tp_pf_pbs; multiconductor=true, solution_builder = get_pbs_solution, kwargs...)
+end
+
+
+""
+function run_tp_pf_pbs(file::String, model_constructor, solver; kwargs...)
+    data = ThreePhasePowerModels.parse_file(file)
+    return PMs.run_generic_model(data, model_constructor, solver, post_tp_pf_pbs; multiconductor=true, solution_builder = get_pbs_solution, kwargs...)
+end
+
+""
+function post_tp_pf_pbs(pm::GenericPowerModel)
+    for c in PMs.conductor_ids(pm)
+        variable_tp_voltage(pm, bounded=false, cnd=c)
+        variable_bus_power_slack(pm, cnd=c)
+        PMs.variable_generation(pm, bounded=false, cnd=c)
+        PMs.variable_branch_flow(pm, bounded=false, cnd=c)
+        PMs.variable_dcline_flow(pm, bounded=false, cnd=c)
+    end
+
+    for c in PMs.conductor_ids(pm)
+        constraint_tp_voltage(pm, cnd=c)
+
+        for (i,bus) in ref(pm, :ref_buses)
+            @assert bus["bus_type"] == 3
+            constraint_tp_theta_ref(pm, i, cnd=c)
+            PMs.constraint_voltage_magnitude_setpoint(pm, i, cnd=c)
+        end
+
+        for (i,bus) in ref(pm, :bus)
+            constraint_kcl_shunt_slack(pm, i, cnd=c)
+
+            # PV Bus Constraints
+            if length(ref(pm, :bus_gens, i)) > 0 && !(i in ids(pm,:ref_buses))
+                # this assumes inactive generators are filtered out of bus_gens
+                @assert bus["bus_type"] == 2
+
+                PMs.constraint_voltage_magnitude_setpoint(pm, i, cnd=c)
+                for j in ref(pm, :bus_gens, i)
+                    PMs.constraint_active_gen_setpoint(pm, j, cnd=c)
+                end
+            end
+        end
+
+        for i in ids(pm, :branch)
+            constraint_ohms_tp_yt_from(pm, i, cnd=c)
+            constraint_ohms_tp_yt_to(pm, i, cnd=c)
+        end
+
+        for (i,dcline) in ref(pm, :dcline)
+            PMs.constraint_active_dcline_setpoint(pm, i, cnd=c)
+
+            f_bus = ref(pm, :bus)[dcline["f_bus"]]
+            if f_bus["bus_type"] == 1
+                PMs.constraint_voltage_magnitude_setpoint(pm, f_bus["index"], cnd=c)
+            end
+
+            t_bus = ref(pm, :bus)[dcline["t_bus"]]
+            if t_bus["bus_type"] == 1
+                PMs.constraint_voltage_magnitude_setpoint(pm, t_bus["index"], cnd=c)
+            end
+        end
+
+    end
+
+    objective_min_bus_power_slack(pm)
+end
+
+
+""
+function get_pbs_solution(pm::GenericPowerModel, sol::Dict{String,Any})
+    PMs.add_bus_voltage_setpoint(sol, pm)
+    PMs.add_generator_power_setpoint(sol, pm)
+    PMs.add_branch_flow_setpoint(sol, pm)
+    add_bus_slack_setpoint(sol, pm)
+end
+
+
+""
+function add_bus_slack_setpoint(sol, pm::GenericPowerModel)
+    PMs.add_setpoint(sol, pm, "bus", "p_slack", :p_slack)
+    PMs.add_setpoint(sol, pm, "bus", "q_slack", :q_slack)
+end
+
+
diff --git a/test/data/matlab/case5_c_m_a.m b/test/data/matlab/case5_c_m_a.m
index b3d935376..c7e4f7c59 100644
--- a/test/data/matlab/case5_c_m_a.m
+++ b/test/data/matlab/case5_c_m_a.m
@@ -37,7 +37,7 @@
 tppmc.gen = [
 	1	 0.0	  40.0	  -30.0	  30.0	 0.0	  40.0	  -30.0	  30.0	 0.0	  40.0	  -30.0	  30.0	  40.0	  30.0	  40.0	  30.0	  40.0	  30.0	 1;
 	1	 0.0	 170.0	 -127.5	 127.5	 0.0 	 170.0	 -127.5	 127.5	 0.0 	 170.0	 -127.5	 127.5	 170.0	 127.5	 170.0	 127.5	 170.0	 127.5	 1;
-	3	 0.0	 520.0	 -390.	 390.0	 0.0 	 520.0	 -390.0	 390.0	 0.0 	 520.0	 -390.0	 390.0	 325.0	 390.0	 325.0	 390.0	 325.0	 390.0	 1;
+	3	 0.0	 520.0	 -390.0	 390.0	 0.0 	 520.0	 -390.0	 390.0	 0.0 	 520.0	 -390.0	 390.0	 325.0	 390.0	 325.0	 390.0	 325.0	 390.0	 1;
 	4	 0.0	 200.0	 -150.0	 150.0	 0.0	 200.0	 -150.0	 150.0	 0.0	 200.0	 -150.0	 150.0	   0.0	 -10.0	   0.0	 -10.0	   0.0	 -10.0	 1;
 	10	 0.0	 600.0	 -450.0	 450.0	 0.0	 600.0	 -450.0	 450.0	 0.0	 600.0	 -450.0	 450.0	 470.0	-165.0	 470.0	-165.0	 470.0	-165.0	 1;
 ];
diff --git a/test/data/matlab/case5_c_m_b.m b/test/data/matlab/case5_c_m_b.m
new file mode 100644
index 000000000..44fc0215a
--- /dev/null
+++ b/test/data/matlab/case5_c_m_b.m
@@ -0,0 +1,62 @@
+% 3 Coupled meshed networks with off diagonal terms in branch Y
+% tests that opf_pbs model converges on an otherwise infeasible model
+
+function tppmc = case5_c_m_b
+tppmc.version = '1'
+
+tppmc.baseMVA = 100.0;
+tppmc.baseKV = 230.0;
+
+%% bus data
+%	bus_i	type	vmin_1	vmax_1	vmin_2	vmax_2	vmin_3	vmax_3	vm_1	va_1	vm_2	va_2	vm_3	va_3
+tppmc.bus = [
+	1	 2	 0.90	 1.10	 0.90	 1.10	 0.90	 1.10	 1.00000	 2.80377	 1.00000	 2.80377	 1.00000	 2.80377;
+	2	 1	 0.90	 1.10	 0.90	 1.10	 0.90	 1.10	 1.08407	-0.73465	 1.08407	-0.73465	 1.08407	-0.73465;
+	3	 2	 0.90	 1.10	 0.90	 1.10	 0.90	 1.10	 1.00000	-0.55972	 1.00000	-0.55972	 1.00000	-0.55972;
+	4	 1	 0.90	 1.10	 0.90	 1.10	 0.90	 1.10	 1.00000	 0.00000	 1.00000	 0.00000	 1.00000	 0.00000;
+	10	 3	 0.90	 1.10	 0.90	 1.10	 0.90	 1.10	 1.00000	 3.59033	 1.00000	 3.59033	 1.00000	 3.59033;
+];
+
+%% load data
+%	load_bus	pd_1	qd_1	pd_2	qd_2	pd_3	qd_3	status
+tppmc.load = [
+	2	 850.0	 100.0	 300.0	 100.0	 300.0	 100.0	 1;
+	3	 300.0	 100.0	 300.0	 300.0	 300.0	 500.0	 1;
+	4	 400.0	 130.0	1150.0	 130.0	 400.0	 130.0	 1;
+];
+
+%% shunt data
+%	shunt_bus	gs_1	bs_1	gs_2	bs_2	gs_3	bs_3	status
+tppmc.shunt = [
+	 4	 2.0	 50.0	 2.0	 50.0	 2.0	 50.0	 1;
+	10	 1.0	 25.0	 1.0	 25.0	 1.0	 25.0	 1;
+];
+
+%% generator data
+%	gen_bus	pmin_1	pmax_1	qmin_1	qmax_1	pmin_2	pmax_2	qmin_2	qmax_2	pmin_3	pmax_3	qmin_3	qmax_3	pg_1	qg_1	pg_2	qg_2	pg_3	qg_3	status
+tppmc.gen = [
+	1	 0.0	  40.0	  -30.0	  30.0	 0.0	  40.0	  -30.0	  30.0	 0.0	  40.0	  -30.0	  30.0	  40.0	  30.0	  40.0	  30.0	  40.0	  30.0	 1;
+	1	 0.0	 170.0	 -127.5	 127.5	 0.0 	 170.0	 -127.5	 127.5	 0.0 	 170.0	 -127.5	 127.5	 170.0	 127.5	 170.0	 127.5	 170.0	 127.5	 1;
+	3	 0.0	 520.0	 -190.0	 190.0	 0.0 	 520.0	 -190.0	 190.0	 0.0 	 520.0	 -190.0	 190.0	 325.0	 390.0	 325.0	 390.0	 325.0	 390.0	 1;
+	10	 0.0	 600.0	 -450.0	 450.0	 0.0	 600.0	 -450.0	 450.0	 0.0	 600.0	 -450.0	 450.0	 470.0	-165.0	 470.0	-165.0	 470.0	-165.0	 1;
+];
+
+%% generator cost data
+%	2	startup	shutdown	n	c(n-1)	...	c0
+tppmc.gencost = [
+	2	 0.0	 0.0	 3	 0.0	 14.0	 0.0;
+	2	 0.0	 0.0	 3	 0.0	 15.0	 0.0;
+	2	 0.0	 0.0	 3	 0.0	 30.0	 0.0;
+	2	 0.0	 0.0	 3	 0.0	 10.0	 0.0;
+];
+
+%% branch data
+%	f_bus	t_bus	r_11	x_11	r_12	x_12	r_13	x_13	r_22	x_22	r_23	x_23	r_33	x_33    b_1    b_2    b_3    rate_a	rate_b	rate_c	angmin	angmax	status
+tppmc.branch = [
+	1	 2	 0.00281	 0.0281	 0.0001	 0.0002	 0.0001   0.0002	 0.00581	 0.0281	 0.0001   0.0002	 0.00281	 0.0281    0.0    0.0    0.0    400	 400	 400	 -30.0	 30.0	 1;
+	1	 4	 0.00304	 0.0304	 0.0001	 0.0002	 0.0001   0.0002	 0.01004	 0.0304	 0.0001   0.0002	 0.00304	 0.0304    0.0    0.0    0.0    426	 426	 426	 -30.0	 30.0	 1;
+	1	 10	 0.00064	 0.0064	 0.0001  0.0002	 0.0001   0.0002	 0.00164	 0.0064	 0.0001   0.0002	 0.00064	 0.0064    0.0    0.0    0.0    426	 426	 426	 -30.0	 30.0	 1;
+	2	 3	 0.00108	 0.0108	 0.0001	 0.0002	 0.0001   0.0002	 0.00208	 0.0108	 0.0001   0.0002	 0.00108	 0.0108    0.0    0.0    0.0    426	 426	 426	 -30.0	 30.0	 1;
+	4	 3	 0.00297	 0.0297	 0.0001	 0.0002	 0.0001   0.0002	 0.00597	 0.0297	 0.0001   0.0002	 0.00297	 0.0297    0.0    0.0    0.0    426	 426	 426	 -30.0	 30.0	 1;
+	4	 10	 0.00297	 0.0297	 0.0001	 0.0002	 0.0001   0.0002	 0.00497	 0.0297	 0.0001   0.0002	 0.00297	 0.0297    0.0    0.0    0.0    240	 240	 240	 -30.0	 30.0	 1;
+];
diff --git a/test/runtests.jl b/test/runtests.jl
index e5599b3ab..43db8bdf4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -40,4 +40,7 @@ pajarito_solver = PajaritoSolver(mip_solver=cbc_solver, cont_solver=ipopt_solver
     include("tp_pf.jl")
 
     include("tp_opf_bf.jl")
+
+    include("tp_debug.jl")
+
 end
diff --git a/test/tp_debug.jl b/test/tp_debug.jl
new file mode 100644
index 000000000..4cd5d670a
--- /dev/null
+++ b/test/tp_debug.jl
@@ -0,0 +1,51 @@
+
+@testset "test ac opf pbs" begin
+
+    @testset "5-bus coupled meshed infeasible case" begin
+        @testset "ac case" begin
+            mp_data = ThreePhasePowerModels.parse_file("../test/data/matlab/case5_c_m_b.m")
+            result = ThreePhasePowerModels.run_tp_opf_pbs(mp_data, PMs.ACPPowerModel, ipopt_solver)
+
+            @test result["status"] == :LocalOptimal
+            @test isapprox(result["objective"], 26.4471; atol = 1e-2)
+        end
+    end
+
+end
+
+@testset "test soc opf pbs" begin
+
+    @testset "5-bus coupled meshed infeasible case" begin
+        mp_data = ThreePhasePowerModels.parse_file("../test/data/matlab/case5_c_m_b.m")
+        result = ThreePhasePowerModels.run_tp_pf_pbs(mp_data, PMs.SOCWRPowerModel, ipopt_solver)
+
+        @test result["status"] == :LocalOptimal
+        @test isapprox(result["objective"], 0.0; atol = 1e-4)
+    end
+
+end
+
+
+@testset "test soc pf pbs" begin
+
+    @testset "5-bus coupled meshed network (b)" begin
+        mp_data = ThreePhasePowerModels.parse_file("../test/data/matlab/case5_c_m_b.m")
+        result = ThreePhasePowerModels.run_tp_pf_pbs(mp_data, PMs.ACPPowerModel, ipopt_solver)
+
+        @test result["status"] == :LocalOptimal
+        @test isapprox(result["objective"], 0.0; atol=1e-4)
+    end
+
+end
+
+@testset "test soc pf pbs" begin
+
+    @testset "5-bus coupled meshed network (b)" begin
+        mp_data = ThreePhasePowerModels.parse_file("../test/data/matlab/case5_c_m_b.m")
+        result = ThreePhasePowerModels.run_tp_pf_pbs(mp_data, PMs.SOCWRPowerModel, ipopt_solver)
+
+        @test result["status"] == :LocalOptimal
+        @test isapprox(result["objective"], 0.0; atol=1e-4)
+    end
+
+end
\ No newline at end of file
diff --git a/test/tp_opf.jl b/test/tp_opf.jl
index ccebc4fec..83def5f50 100644
--- a/test/tp_opf.jl
+++ b/test/tp_opf.jl
@@ -134,6 +134,23 @@ end
             @test isapprox(result["solution"]["gen"]["1"]["qg"][3], 0.0357746; atol = 1e-3)
         end
     end
+    @testset "5-bus coupled meshed infeasible case" begin
+        @testset "ac case" begin
+            mp_data = ThreePhasePowerModels.parse_file("../test/data/matlab/case5_c_m_b.m")
+            result = run_tp_opf(mp_data, PMs.ACPPowerModel, ipopt_solver)
+
+            @test result["status"] == :LocalInfeasible
+        end
+        #=
+        # omit due to large number of terminal warnings
+        @testset "soc case" begin
+            mp_data = ThreePhasePowerModels.parse_file("../test/data/matlab/case5_c_m_b.m")
+            result = run_tp_opf(mp_data, PMs.SOCWRPowerModel, ipopt_solver)
+
+            @test result["status"] == :LocalInfeasible
+        end
+        =#
+    end
     @testset "5-bus coupled radial no shunt case" begin
         mp_data = ThreePhasePowerModels.parse_file("../test/data/matlab/case5_c_r_a.m")
         result = run_tp_opf(mp_data, PMs.ACPPowerModel, ipopt_solver)
diff --git a/test/tp_pf.jl b/test/tp_pf.jl
index 67b83e90f..610a51ee1 100644
--- a/test/tp_pf.jl
+++ b/test/tp_pf.jl
@@ -6,13 +6,20 @@
         @test isapprox(result["objective"], 0.0; atol=1e-2)
     end
 
-    @testset "5-bus coupled meshed network" begin
+    @testset "5-bus coupled meshed network (a)" begin
         result = run_ac_tp_pf("../test/data/matlab/case5_c_m_a.m", ipopt_solver)
 
         @test result["status"] == :LocalOptimal
         @test isapprox(result["objective"], 0.0; atol=1e-2)
     end
 
+    @testset "5-bus coupled meshed network (b)" begin
+        result = run_ac_tp_pf("../test/data/matlab/case5_c_m_b.m", ipopt_solver)
+
+        @test result["status"] == :LocalOptimal
+        @test isapprox(result["objective"], 0.0; atol=1e-2)
+    end
+
     @testset "5-bus independent radial w/ shunts" begin
         result = run_ac_tp_pf("../test/data/matlab/case5_i_r_b.m", ipopt_solver)
 
@@ -65,3 +72,13 @@
         @test isapprox(sum(sol["solution"]["gen"]["1"]["qg"] * sol["solution"]["baseMVA"]), 0.00932693; atol=1e-4)
     end
 end
+
+
+@testset "test soc pf" begin
+    @testset "5-bus coupled meshed network (b)" begin
+        result = run_tp_pf("../test/data/matlab/case5_c_m_b.m", PMs.SOCWRPowerModel, ipopt_solver)
+
+        @test result["status"] == :LocalOptimal
+        @test isapprox(result["objective"], 0.0; atol=1e-2)
+    end
+end