diff --git a/src/form_decomposition/boundary.jl b/src/form_decomposition/boundary.jl
index 52144ae..1e806f0 100644
--- a/src/form_decomposition/boundary.jl
+++ b/src/form_decomposition/boundary.jl
@@ -117,6 +117,25 @@ function constraint_boundary_voltage_magnitude(pmd::_PMD.AbstractUnbalancedPower
 end
 
 
+"""
+    function constraint_boundary_voltage_angle(
+        pmd::_PMD.AbstractUnbalancedPowerModel,
+        i::Int;
+        nw::Int=nw_id_default
+    )
+
+General voltage angle boundary constraint.
+"""
+function constraint_boundary_voltage_angle(pmd::_PMD.AbstractUnbalancedPowerModel, i::Int; nw::Int=nw_id_default)
+
+    boundary = _PMD.ref(pmd, nw, :boundary, i)
+    t_bus = boundary["t_bus"] # convention: to bus Distribution always!
+
+    constraint_boundary_voltage_angle(pmd, i, t_bus, boundary["f_connections"], boundary["t_connections"]; nw=nw)
+
+end
+
+
 """
     function constraint_transmission_power_balance_boundary(
         pm::_PM.AbstractPowerModel,
diff --git a/src/form_decomposition/ivr.jl b/src/form_decomposition/ivr.jl
index 50820c6..01b8e0e 100644
--- a/src/form_decomposition/ivr.jl
+++ b/src/form_decomposition/ivr.jl
@@ -53,7 +53,7 @@ end
 
 """
     function constraint_boundary_voltage_magnitude(
-        pmd::_PMD.ACPUPowerModel,
+        pmd::_PMD.IVRUPowerModel,
         i::Int,
         f_idx::Tuple{Int,Int,Int},
         f_connections::Vector{Int},
@@ -75,6 +75,35 @@ function constraint_boundary_voltage_magnitude(pmd::_PMD.IVRUPowerModel, ::Int,
 end
 
 
+"""
+    function constraint_boundary_voltage_angle(
+        pmd::_PMD.IVRUPowerModel,
+        i::Int,
+        f_idx::Tuple{Int,Int,Int},
+        f_connections::Vector{Int},
+        t_connections::Vector{Int};
+        nw::Int = nw_id_default
+    )
+
+IVRU boundary bus voltage angle constraints.
+"""
+function constraint_boundary_voltage_angle(pmd::_PMD.IVRUPowerModel, ::Int, t_bus::Int, ::Vector{Int}, ::Vector{Int}; nw::Int=nw_id_default)
+
+    vr_source = _PMD.var(pmd, nw, :vr, t_bus)
+    vi_source = _PMD.var(pmd, nw, :vi, t_bus)
+
+    # Add constraint(s): r->real, i->imaginary
+    JuMP.@constraint(pmd.model, vi_source[1] == 0.0)
+
+    # Add constraints related to 120 degrees offset for the distribution b and c phases
+    shift_120degs_rad = deg2rad(120)
+
+    # TODO: These are non-linear constraints due to transformation to degrees of phase a angle (another way - non-linear may be possible)
+    JuMP.@NLconstraint(pmd.model, vi_source[2] == tan(atan(vi_source[1]/vr_source[1]) - shift_120degs_rad)*vr_source[2])
+    JuMP.@NLconstraint(pmd.model, vi_source[3] == tan(atan(vi_source[1]/vr_source[1]) + shift_120degs_rad)*vr_source[3])
+
+end
+
 
 """
     function constraint_transmission_boundary_power_shared_vars_scaled(
diff --git a/src/prob/opfitd_decomposition.jl b/src/prob/opfitd_decomposition.jl
index 7945cd7..6513103 100755
--- a/src/prob/opfitd_decomposition.jl
+++ b/src/prob/opfitd_decomposition.jl
@@ -476,6 +476,7 @@ function build_opfitd_decomposition(pmd_model::_PMD.AbstractUnbalancedIVRModel)
     for i in _PMD.ids(pmd_model, :boundary)
         constraint_boundary_power(pmd_model, i)
         constraint_boundary_voltage_magnitude(pmd_model, i)
+        constraint_boundary_voltage_angle(pmd_model, i)
     end
 
     # PMD objective