diff --git a/CHANGELOG.md b/CHANGELOG.md
index 67f0d54be..9279c5b9c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,7 +2,7 @@
 
 ## staged
 
-- none
+- Added UBF matrix power variables for switches [#423](https://github.com/lanl-ansi/PowerModelsDistribution.jl/issues/423)
 
 ## v0.14.5
 
diff --git a/src/form/bf_mx.jl b/src/form/bf_mx.jl
index a69886d5a..524b5bef3 100644
--- a/src/form/bf_mx.jl
+++ b/src/form/bf_mx.jl
@@ -137,6 +137,70 @@ function variable_mc_branch_power(pm::AbstractUBFModels; nw::Int=nw_id_default,
 end
 
 
+"matrix power variables for switches"
+function variable_mc_switch_power(pm::AbstractUBFModels; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
+    # calculate S bound
+    switch_arcs = Vector{Tuple{Int,Int,Int}}(ref(pm, nw, :arcs_switch))
+    connections = Dict{Tuple{Int,Int,Int},Vector{Int}}((l,i,j) => connections for (bus,entry) in ref(pm, nw, :bus_arcs_conns_switch) for ((l,i,j), connections) in entry)
+
+    if bounded
+        bound = Dict{eltype(switch_arcs), Matrix{Real}}()
+        for (l, switch) in ref(pm, nw, :switch)
+            bus_fr = ref(pm, nw, :bus, switch["f_bus"])
+            bus_to = ref(pm, nw, :bus, switch["t_bus"])
+
+            smax_fr = _calc_branch_power_max(switch, bus_fr)
+            smax_to = _calc_branch_power_max(switch, bus_to)
+            cmax_fr, cmax_to = _calc_branch_current_max_frto(switch, bus_fr, bus_to)
+
+            tuple_fr = (l, bus_fr["index"], bus_to["index"])
+            tuple_to = (l, bus_to["index"], bus_fr["index"])
+
+            bound[tuple_fr] = bus_fr["vmax"][[findfirst(isequal(c), bus_fr["terminals"]) for c in switch["f_connections"]]].*cmax_fr'
+            bound[tuple_to] = bus_to["vmax"][[findfirst(isequal(c), bus_to["terminals"]) for c in switch["t_connections"]]].*cmax_to'
+
+            for (idx, (fc,tc)) in enumerate(zip(switch["f_connections"], switch["t_connections"]))
+                bound[tuple_fr][idx,idx] = smax_fr[idx]
+                bound[tuple_to][idx,idx] = smax_to[idx]
+            end
+        end
+        # create matrix variables
+        (P,Q) = variable_mx_complex(pm.model, switch_arcs, connections, connections; symm_bound=bound, name=("Psw", "Qsw"), prefix="$nw")
+    else
+        (P,Q) = variable_mx_complex(pm.model, switch_arcs, connections, connections; name=("Psw", "Qsw"), prefix="$nw")
+    end
+
+    # this explicit type erasure is necessary
+    P_expr = merge(
+        Dict{Any,Any}( (l,i,j) => P[(l,i,j)] for (l,i,j) in ref(pm, nw, :arcs_switch_from) ),
+        Dict( (l,j,i) => -1.0.*P[(l,i,j)] for (l,i,j) in ref(pm, nw, :arcs_switch_from))
+    )
+    Q_expr = merge(
+        Dict{Any,Any}( (l,i,j) => Q[(l,i,j)] for (l,i,j) in ref(pm, nw, :arcs_switch_from) ),
+        Dict( (l,j,i) => -1.0*Q[(l,i,j)] for (l,i,j) in ref(pm, nw, :arcs_switch_from))
+    )
+
+    # This is needed to get around error: "unexpected affine expression in nlconstraint"
+    (P_aux,Q_aux) = variable_mx_complex(pm.model, switch_arcs, connections, connections; name=("Psw_aux", "Qsw_aux"), prefix="$nw")
+    for (l,i,j) in switch_arcs
+        JuMP.@constraint(pm.model, P_expr[(l,i,j)] .== P_aux[(l,i,j)])
+        JuMP.@constraint(pm.model, Q_expr[(l,i,j)] .== Q_aux[(l,i,j)])
+    end
+
+    # save reference
+    var(pm, nw)[:Psw] = P_aux
+    var(pm, nw)[:Qsw] = Q_aux
+
+    var(pm, nw)[:psw] = Dict([(id,LinearAlgebra.diag(P_aux[id])) for id in switch_arcs])
+    var(pm, nw)[:qsw] = Dict([(id,LinearAlgebra.diag(Q_aux[id])) for id in switch_arcs])
+
+    report && _IM.sol_component_value_edge(pm, pmd_it_sym, nw, :switch, :Pf, :Pt, ref(pm, nw, :arcs_switch_from), ref(pm, nw, :arcs_switch_to), P_expr)
+    report && _IM.sol_component_value_edge(pm, pmd_it_sym, nw, :switch, :Qf, :Qt, ref(pm, nw, :arcs_switch_from), ref(pm, nw, :arcs_switch_to), Q_expr)
+    report && _IM.sol_component_value_edge(pm, pmd_it_sym, nw, :switch, :pf, :pt, ref(pm, nw, :arcs_switch_from), ref(pm, nw, :arcs_switch_to), Dict([(id,LinearAlgebra.diag(P_expr[id])) for id in switch_arcs]))
+    report && _IM.sol_component_value_edge(pm, pmd_it_sym, nw, :switch, :qf, :qt, ref(pm, nw, :arcs_switch_from), ref(pm, nw, :arcs_switch_to), Dict([(id,LinearAlgebra.diag(Q_expr[id])) for id in switch_arcs]))
+end
+
+
 "defines matrix transformer power variables for the unbalanced branch flow models"
 function variable_mc_transformer_power(pm::AbstractUBFModels; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
     transformer_arcs = Vector{Tuple{Int,Int,Int}}(ref(pm, nw, :arcs_transformer))
diff --git a/test/opf_bf.jl b/test/opf_bf.jl
index f3eb4282d..ced221393 100644
--- a/test/opf_bf.jl
+++ b/test/opf_bf.jl
@@ -169,5 +169,56 @@
                 @test isapprox(result["objective"], 21.17; atol = 5e-2)
             end
         end
+
+        @testset "test ubf relaxations with with switches" begin
+            data = deepcopy(case3_balanced_switch)
+            make_lossless!(data; exclude=["line", "linecode"])
+            remove_line_limits!(data)
+            apply_voltage_bounds!(data; vm_lb=0.9, vm_ub=1.1)
+
+            data["settings"]["sbase_default"] = 1.0
+
+            data["generator"] = Dict{String,Any}(
+                "1" => Dict{String,Any}(
+                    "bus" => "primary",
+                    "connections" => [1, 2, 3, 4],
+                    "cost_pg_parameters" => [0.0, 1200.0, 0.0],
+                    "qg_lb" => fill(0.0, 3),
+                    "qg_ub" => fill(0.0, 3),
+                    "pg_ub" => fill(10, 3),
+                    "pg_lb" => fill(0, 3),
+                    "configuration" => WYE,
+                    "status" => ENABLED
+                )
+            )
+
+            merge!(data["voltage_source"]["source"], Dict{String,Any}(
+                "cost_pg_parameters" => [0.0, 1000.0, 0.0],
+                "pg_lb" => fill(  0.0, 3),
+                "pg_ub" => fill( 10.0, 3),
+                "qg_lb" => fill(-10.0, 3),
+                "qg_ub" => fill( 10.0, 3),
+            ))
+
+            for (_,line) in data["line"]
+                line["sm_ub"] = fill(10.0, 3)
+            end
+
+            @testset "test SOCNLPUBF opf with switches" begin
+                result = solve_mc_opf(data, SOCNLPUBFPowerModel, ipopt_solver)
+
+                @test result["termination_status"] == LOCALLY_SOLVED || result["termination_status"] == ALMOST_LOCALLY_SOLVED
+                @test all(isapprox.(result["solution"]["switch"]["ohline"]["pf"], [6.0, 6.0, 6.0]; atol=1e-1))
+                @test isapprox(result["objective"], 18.1824; atol=2e-1)
+            end
+
+            @testset "test SOCConicUBF opf with switches" begin
+                result = solve_mc_opf(data, SOCConicUBFPowerModel, scs_solver)
+
+                @test result["termination_status"] == OPTIMAL || result["termination_status"] == ALMOST_OPTIMAL
+                @test all(isapprox.(result["solution"]["switch"]["ohline"]["pf"], [6.0, 6.0, 6.0]; atol=1e-1))
+                @test isapprox(result["objective"], 18.1824; atol=2e-1)
+            end
+        end
     end
 end
diff --git a/test/test_cases.jl b/test/test_cases.jl
index 43bd2dee9..925b43f96 100644
--- a/test/test_cases.jl
+++ b/test/test_cases.jl
@@ -12,6 +12,7 @@ case3_balanced_pv = parse_file("../test/data/opendss/case3_balanced_pv.dss")
 case3_unbalanced_1phase_pv = parse_file("../test/data/opendss/case3_unbalanced_1phase-pv.dss")
 case3_blanced_basefreq = parse_file("../test/data/opendss/case3_balanced_basefreq.dss")
 case3_balanced_battery = parse_file("../test/data/opendss/case3_balanced_battery.dss")
+case3_balanced_switch = parse_file("../test/data/opendss/case3_balanced_switch.dss")
 
 # 3-bus unbalanced
 case3_unbalanced = parse_file("../test/data/opendss/case3_unbalanced.dss")