diff --git a/benchmark/run_benchmarks.jl b/benchmark/run_benchmarks.jl
index fd4b2943..bdc66680 100644
--- a/benchmark/run_benchmarks.jl
+++ b/benchmark/run_benchmarks.jl
@@ -37,7 +37,7 @@ end
 
 function poi_add_parameters_and_variables_alternating(N::Int)
     model = POI.Optimizer(MOI.Utilities.Model{Float64}())
-    for i in 1:Int(N / 2)
+    for i in 1:Int(N/2)
         MOI.add_variable(model)
         MOI.add_constrained_variable(model, MOI.Parameter(1.0))
     end
@@ -73,13 +73,9 @@ end
 function poi_add_saf_variables_and_parameters_ctr(N::Int, M::Int)
     model = POI.Optimizer(MOI.Utilities.Model{Float64}())
     x = MOI.add_variables(model, N / 2)
-    y =
-        first.(
-            MOI.add_constrained_variable.(
-                model,
-                MOI.Parameter.(ones(Int(N / 2))),
-            )
-        )
+    y = first.(
+        MOI.add_constrained_variable.(model, MOI.Parameter.(ones(Int(N / 2)))),
+    )
     for _ in 1:M
         MOI.add_constraint(
             model,
@@ -99,13 +95,9 @@ function poi_add_saf_variables_and_parameters_ctr_parameter_update(
 )
     model = POI.Optimizer(MOI.Utilities.Model{Float64}())
     x = MOI.add_variables(model, N / 2)
-    y =
-        first.(
-            MOI.add_constrained_variable.(
-                model,
-                MOI.Parameter.(ones(Int(N / 2))),
-            )
-        )
+    y = first.(
+        MOI.add_constrained_variable.(model, MOI.Parameter.(ones(Int(N / 2)))),
+    )
     for _ in 1:M
         MOI.add_constraint(
             model,
@@ -158,13 +150,9 @@ end
 function poi_add_sqf_variables_parameters_ctr(N::Int, M::Int)
     model = POI.Optimizer(MOI.Utilities.Model{Float64}())
     x = MOI.add_variables(model, N / 2)
-    y =
-        first.(
-            MOI.add_constrained_variable.(
-                model,
-                MOI.Parameter.(ones(Int(N / 2))),
-            )
-        )
+    y = first.(
+        MOI.add_constrained_variable.(model, MOI.Parameter.(ones(Int(N / 2)))),
+    )
     for _ in 1:M
         MOI.add_constraint(
             model,
@@ -182,13 +170,9 @@ end
 function poi_add_sqf_variables_parameters_ctr_parameter_update(N::Int, M::Int)
     model = POI.Optimizer(MOI.Utilities.Model{Float64}())
     x = MOI.add_variables(model, N / 2)
-    y =
-        first.(
-            MOI.add_constrained_variable.(
-                model,
-                MOI.Parameter.(ones(Int(N / 2))),
-            )
-        )
+    y = first.(
+        MOI.add_constrained_variable.(model, MOI.Parameter.(ones(Int(N / 2)))),
+    )
     for _ in 1:M
         MOI.add_constraint(
             model,
@@ -208,13 +192,9 @@ end
 function poi_add_sqf_parameters_parameters_ctr(N::Int, M::Int)
     model = POI.Optimizer(MOI.Utilities.Model{Float64}())
     x = MOI.add_variables(model, N / 2)
-    y =
-        first.(
-            MOI.add_constrained_variable.(
-                model,
-                MOI.Parameter.(ones(Int(N / 2))),
-            )
-        )
+    y = first.(
+        MOI.add_constrained_variable.(model, MOI.Parameter.(ones(Int(N / 2)))),
+    )
     for _ in 1:M
         MOI.add_constraint(
             model,
@@ -232,13 +212,9 @@ end
 function poi_add_sqf_parameters_parameters_ctr_parameter_update(N::Int, M::Int)
     model = POI.Optimizer(MOI.Utilities.Model{Float64}())
     x = MOI.add_variables(model, N / 2)
-    y =
-        first.(
-            MOI.add_constrained_variable.(
-                model,
-                MOI.Parameter.(ones(Int(N / 2))),
-            )
-        )
+    y = first.(
+        MOI.add_constrained_variable.(model, MOI.Parameter.(ones(Int(N / 2)))),
+    )
     for _ in 1:M
         MOI.add_constraint(
             model,
@@ -284,13 +260,9 @@ end
 function poi_add_saf_variables_and_parameters_obj(N::Int, M::Int)
     model = POI.Optimizer(MOI.Utilities.Model{Float64}())
     x = MOI.add_variables(model, N / 2)
-    y =
-        first.(
-            MOI.add_constrained_variable.(
-                model,
-                MOI.Parameter.(ones(Int(N / 2))),
-            )
-        )
+    y = first.(
+        MOI.add_constrained_variable.(model, MOI.Parameter.(ones(Int(N / 2)))),
+    )
     for _ in 1:M
         MOI.set(
             model,
@@ -310,13 +282,9 @@ function poi_add_saf_variables_and_parameters_obj_parameter_update(
 )
     model = POI.Optimizer(MOI.Utilities.Model{Float64}())
     x = MOI.add_variables(model, N / 2)
-    y =
-        first.(
-            MOI.add_constrained_variable.(
-                model,
-                MOI.Parameter.(ones(Int(N / 2))),
-            )
-        )
+    y = first.(
+        MOI.add_constrained_variable.(model, MOI.Parameter.(ones(Int(N / 2)))),
+    )
     for _ in 1:M
         MOI.set(
             model,
@@ -371,13 +339,9 @@ end
 function poi_add_sqf_variables_parameters_obj(N::Int, M::Int)
     model = POI.Optimizer(MOI.Utilities.Model{Float64}())
     x = MOI.add_variables(model, N / 2)
-    y =
-        first.(
-            MOI.add_constrained_variable.(
-                model,
-                MOI.Parameter.(ones(Int(N / 2))),
-            )
-        )
+    y = first.(
+        MOI.add_constrained_variable.(model, MOI.Parameter.(ones(Int(N / 2)))),
+    )
     for _ in 1:M
         MOI.set(
             model,
@@ -395,13 +359,9 @@ end
 function poi_add_sqf_variables_parameters_obj_parameter_update(N::Int, M::Int)
     model = POI.Optimizer(MOI.Utilities.Model{Float64}())
     x = MOI.add_variables(model, N / 2)
-    y =
-        first.(
-            MOI.add_constrained_variable.(
-                model,
-                MOI.Parameter.(ones(Int(N / 2))),
-            )
-        )
+    y = first.(
+        MOI.add_constrained_variable.(model, MOI.Parameter.(ones(Int(N / 2)))),
+    )
     for _ in 1:M
         MOI.set(
             model,
@@ -423,13 +383,9 @@ end
 function poi_add_sqf_parameters_parameters_obj(N::Int, M::Int)
     model = POI.Optimizer(MOI.Utilities.Model{Float64}())
     x = MOI.add_variables(model, N / 2)
-    y =
-        first.(
-            MOI.add_constrained_variable.(
-                model,
-                MOI.Parameter.(ones(Int(N / 2))),
-            )
-        )
+    y = first.(
+        MOI.add_constrained_variable.(model, MOI.Parameter.(ones(Int(N / 2)))),
+    )
     for _ in 1:M
         MOI.set(
             model,
@@ -447,13 +403,9 @@ end
 function poi_add_sqf_parameters_parameters_obj_parameter_update(N::Int, M::Int)
     model = POI.Optimizer(MOI.Utilities.Model{Float64}())
     x = MOI.add_variables(model, N / 2)
-    y =
-        first.(
-            MOI.add_constrained_variable.(
-                model,
-                MOI.Parameter.(ones(Int(N / 2))),
-            )
-        )
+    y = first.(
+        MOI.add_constrained_variable.(model, MOI.Parameter.(ones(Int(N / 2)))),
+    )
     for _ in 1:M
         MOI.set(
             model,
diff --git a/benchmark/run_benchmarks_jump.jl b/benchmark/run_benchmarks_jump.jl
index fe95b743..877b2707 100644
--- a/benchmark/run_benchmarks_jump.jl
+++ b/benchmark/run_benchmarks_jump.jl
@@ -17,7 +17,7 @@ function moi_add_variables(N::Int)
             MOI.Utilities.AUTOMATIC,
         ),
     )
-    @variable(model, x[i = 1:N])
+    @variable(model, x[i=1:N])
     return nothing
 end
 
@@ -30,7 +30,7 @@ function poi_add_variables(N::Int)
             ),
         ),
     )
-    @variable(model, x[i = 1:N])
+    @variable(model, x[i=1:N])
     return nothing
 end
 
@@ -43,7 +43,7 @@ function poi_add_parameters(N::Int)
             ),
         ),
     )
-    @variable(model, x[i = 1:N] in MOI.Parameter(1.0))
+    @variable(model, x[i=1:N] in MOI.Parameter(1.0))
     return nothing
 end
 
@@ -56,7 +56,7 @@ function poi_add_parameters_and_variables(N::Int)
             ),
         ),
     )
-    @variable(model, x[i = 1:N] in MOI.Parameter(1.0))
+    @variable(model, x[i=1:N] in MOI.Parameter(1.0))
     return nothing
 end
 
@@ -69,7 +69,7 @@ function poi_add_parameters_and_variables_alternating(N::Int)
             ),
         ),
     )
-    for i in 1:Int(N / 2)
+    for i in 1:Int(N/2)
         @variable(model)
         @variable(model, set = MOI.Parameter(1.0))
     end
@@ -83,8 +83,8 @@ function moi_add_saf_ctr(N::Int, M::Int)
             MOI.Utilities.AUTOMATIC,
         ),
     )
-    @variable(model, x[i = 1:N])
-    @constraint(model, cons[i = 1:M], sum(x) >= 1)
+    @variable(model, x[i=1:N])
+    @constraint(model, cons[i=1:M], sum(x) >= 1)
     return nothing
 end
 
@@ -97,8 +97,8 @@ function poi_add_saf_ctr(N::Int, M::Int)
             ),
         ),
     )
-    @variable(model, x[i = 1:N])
-    @constraint(model, cons[i = 1:M], sum(x) >= 1)
+    @variable(model, x[i=1:N])
+    @constraint(model, cons[i=1:M], sum(x) >= 1)
     return nothing
 end
 
@@ -111,9 +111,9 @@ function poi_add_saf_variables_and_parameters_ctr(N::Int, M::Int)
             ),
         ),
     )
-    @variable(model, x[i = 1:Int(N / 2)])
-    @variable(model, p[i = 1:Int(N / 2)] in MOI.Parameter.(0))
-    @constraint(model, con[i = 1:M], sum(x) + sum(p) >= 1)
+    @variable(model, x[i=1:Int(N/2)])
+    @variable(model, p[i=1:Int(N/2)] in MOI.Parameter.(0))
+    @constraint(model, con[i=1:M], sum(x) + sum(p) >= 1)
     return nothing
 end
 
@@ -129,9 +129,9 @@ function poi_add_saf_variables_and_parameters_ctr_parameter_update(
             ),
         ),
     )
-    @variable(model, x[i = 1:Int(N / 2)])
-    @variable(model, p[i = 1:Int(N / 2)] in MOI.Parameter.(0))
-    @constraint(model, con[i = 1:M], sum(x) + sum(p) >= 1)
+    @variable(model, x[i=1:Int(N/2)])
+    @variable(model, p[i=1:Int(N/2)] in MOI.Parameter.(0))
+    @constraint(model, con[i=1:M], sum(x) + sum(p) >= 1)
     MOI.set.(model, POI.ParameterValue(), p, 0.5)
     POI.update_parameters!(backend(model))
     return nothing
@@ -144,8 +144,8 @@ function moi_add_sqf_variables_ctr(N::Int, M::Int)
             MOI.Utilities.AUTOMATIC,
         ),
     )
-    @variable(model, x[i = 1:N])
-    @constraint(model, con[i = 1:M], dot(x, x) >= 1)
+    @variable(model, x[i=1:N])
+    @constraint(model, con[i=1:M], dot(x, x) >= 1)
     return nothing
 end
 
@@ -158,8 +158,8 @@ function poi_add_sqf_variables_ctr(N::Int, M::Int)
             ),
         ),
     )
-    @variable(model, x[i = 1:N])
-    @constraint(model, con[i = 1:M], dot(x, x) >= 1)
+    @variable(model, x[i=1:N])
+    @constraint(model, con[i=1:M], dot(x, x) >= 1)
     return nothing
 end
 
@@ -172,9 +172,9 @@ function poi_add_sqf_variables_parameters_ctr(N::Int, M::Int)
             ),
         ),
     )
-    @variable(model, x[i = 1:Int(N / 2)])
-    @variable(model, p[i = 1:Int(N / 2)] in MOI.Parameter.(1))
-    @constraint(model, con[i = 1:M], dot(x, p) >= 1)
+    @variable(model, x[i=1:Int(N/2)])
+    @variable(model, p[i=1:Int(N/2)] in MOI.Parameter.(1))
+    @constraint(model, con[i=1:M], dot(x, p) >= 1)
     return nothing
 end
 
@@ -187,9 +187,9 @@ function poi_add_sqf_variables_parameters_ctr_parameter_update(N::Int, M::Int)
             ),
         ),
     )
-    @variable(model, x[i = 1:Int(N / 2)])
-    @variable(model, p[i = 1:Int(N / 2)] in MOI.Parameter.(1))
-    @constraint(model, con[i = 1:M], dot(x, p) >= 1)
+    @variable(model, x[i=1:Int(N/2)])
+    @variable(model, p[i=1:Int(N/2)] in MOI.Parameter.(1))
+    @constraint(model, con[i=1:M], dot(x, p) >= 1)
     MOI.set.(model, POI.ParameterValue(), p, 0.5)
     POI.update_parameters!(backend(model))
     return nothing
@@ -204,9 +204,9 @@ function poi_add_sqf_parameters_parameters_ctr(N::Int, M::Int)
             ),
         ),
     )
-    @variable(model, x[i = 1:Int(N / 2)])
-    @variable(model, p[i = 1:Int(N / 2)] in MOI.Parameter.(1))
-    @constraint(model, con[i = 1:M], dot(p, p) >= 1)
+    @variable(model, x[i=1:Int(N/2)])
+    @variable(model, p[i=1:Int(N/2)] in MOI.Parameter.(1))
+    @constraint(model, con[i=1:M], dot(p, p) >= 1)
     return nothing
 end
 
@@ -219,9 +219,9 @@ function poi_add_sqf_parameters_parameters_ctr_parameter_update(N::Int, M::Int)
             ),
         ),
     )
-    @variable(model, x[i = 1:Int(N / 2)])
-    @variable(model, p[i = 1:Int(N / 2)] in MOI.Parameter.(1))
-    @constraint(model, con[i = 1:M], dot(p, p) >= 1)
+    @variable(model, x[i=1:Int(N/2)])
+    @variable(model, p[i=1:Int(N/2)] in MOI.Parameter.(1))
+    @constraint(model, con[i=1:M], dot(p, p) >= 1)
     MOI.set.(model, POI.ParameterValue(), p, 0.5)
     POI.update_parameters!(backend(model))
     return nothing
@@ -234,7 +234,7 @@ function moi_add_saf_obj(N::Int, M::Int)
             MOI.Utilities.AUTOMATIC,
         ),
     )
-    @variable(model, x[i = 1:N])
+    @variable(model, x[i=1:N])
     @objective(model, Min, sum(x))
     return nothing
 end
@@ -248,7 +248,7 @@ function poi_add_saf_obj(N::Int, M::Int)
             ),
         ),
     )
-    @variable(model, x[i = 1:N])
+    @variable(model, x[i=1:N])
     for _ in 1:M
         @objective(model, Min, sum(x))
     end
@@ -264,8 +264,8 @@ function poi_add_saf_variables_and_parameters_obj(N::Int, M::Int)
             ),
         ),
     )
-    @variable(model, x[i = 1:Int(N / 2)])
-    @variable(model, p[i = 1:Int(N / 2)] in MOI.Parameter.(1))
+    @variable(model, x[i=1:Int(N/2)])
+    @variable(model, p[i=1:Int(N/2)] in MOI.Parameter.(1))
     for _ in 1:M
         @objective(model, Min, sum(x) + sum(p))
     end
@@ -284,8 +284,8 @@ function poi_add_saf_variables_and_parameters_obj_parameter_update(
             ),
         ),
     )
-    @variable(model, x[i = 1:Int(N / 2)])
-    @variable(model, p[i = 1:Int(N / 2)] in MOI.Parameter.(1))
+    @variable(model, x[i=1:Int(N/2)])
+    @variable(model, p[i=1:Int(N/2)] in MOI.Parameter.(1))
     for _ in 1:M
         @objective(model, Min, sum(x) + sum(p))
     end
@@ -303,7 +303,7 @@ function moi_add_sqf_variables_obj(N::Int, M::Int)
             MOI.Utilities.AUTOMATIC,
         ),
     )
-    @variable(model, x[i = 1:N])
+    @variable(model, x[i=1:N])
     for _ in 1:M
         @objective(model, Min, dot(x, x))
     end
@@ -319,7 +319,7 @@ function poi_add_sqf_variables_obj(N::Int, M::Int)
             ),
         ),
     )
-    @variable(model, x[i = 1:N])
+    @variable(model, x[i=1:N])
     for _ in 1:M
         @objective(model, Min, dot(x, x))
     end
@@ -335,8 +335,8 @@ function poi_add_sqf_variables_parameters_obj(N::Int, M::Int)
             ),
         ),
     )
-    @variable(model, x[i = 1:Int(N / 2)])
-    @variable(model, p[i = 1:Int(N / 2)] in MOI.Parameter.(1))
+    @variable(model, x[i=1:Int(N/2)])
+    @variable(model, p[i=1:Int(N/2)] in MOI.Parameter.(1))
     for _ in 1:M
         @objective(model, Min, dot(x, p))
     end
@@ -352,8 +352,8 @@ function poi_add_sqf_variables_parameters_obj_parameter_update(N::Int, M::Int)
             ),
         ),
     )
-    @variable(model, x[i = 1:Int(N / 2)])
-    @variable(model, p[i = 1:Int(N / 2)] in MOI.Parameter.(1))
+    @variable(model, x[i=1:Int(N/2)])
+    @variable(model, p[i=1:Int(N/2)] in MOI.Parameter.(1))
     for _ in 1:M
         @objective(model, Min, dot(x, p))
     end
@@ -373,8 +373,8 @@ function poi_add_sqf_parameters_parameters_obj(N::Int, M::Int)
             ),
         ),
     )
-    @variable(model, x[i = 1:Int(N / 2)])
-    @variable(model, p[i = 1:Int(N / 2)] in MOI.Parameter.(1))
+    @variable(model, x[i=1:Int(N/2)])
+    @variable(model, p[i=1:Int(N/2)] in MOI.Parameter.(1))
     for _ in 1:M
         @objective(model, Min, dot(p, p))
     end
@@ -390,8 +390,8 @@ function poi_add_sqf_parameters_parameters_obj_parameter_update(N::Int, M::Int)
             ),
         ),
     )
-    @variable(model, x[i = 1:Int(N / 2)])
-    @variable(model, p[i = 1:Int(N / 2)] in MOI.Parameter.(1))
+    @variable(model, x[i=1:Int(N/2)])
+    @variable(model, p[i=1:Int(N/2)] in MOI.Parameter.(1))
     for _ in 1:M
         @objective(model, Min, dot(p, p))
     end
diff --git a/benchmark/run_benders_quantile_regression_benchmark.jl b/benchmark/run_benders_quantile_regression_benchmark.jl
index 7765efec..5210908a 100644
--- a/benchmark/run_benders_quantile_regression_benchmark.jl
+++ b/benchmark/run_benders_quantile_regression_benchmark.jl
@@ -139,10 +139,10 @@ function slave_model(PARAM, K)
         # Variables are added to the optimization model, while parameters
         # are not. Parameters are merged with LP problem constants and do not
         # increase the model dimensions.
-        @variable(slave, β[i = 1:N_Candidates] == 0, Param())
+        @variable(slave, β[i=1:N_Candidates] == 0, Param())
     elseif PARAM == 1
         # Create parameters
-        @variable(slave, β[i = 1:N_Candidates] in MOI.Parameter.(0.0))
+        @variable(slave, β[i=1:N_Candidates] in MOI.Parameter.(0.0))
     else
         # Create fixed variables
         @variables(slave, begin
diff --git a/src/parametric_functions.jl b/src/parametric_functions.jl
index d5932358..2a10ed9b 100644
--- a/src/parametric_functions.jl
+++ b/src/parametric_functions.jl
@@ -705,60 +705,92 @@ function _parametric_affine_terms(
             base + term.scalar_term.coefficient * model.parameters[p_idx_val]
     end
 
-    # TODO: check if affine data should only contains variables that appear in pv
-    for (var, coef) in f.affine_data
-        if !haskey(param_terms_dict, var)
-            param_terms_dict[var] = zero(T)
+    # affine data only contain variables that appear in pv
+    for ((var, output_idx), coef) in f.affine_data
+        if !haskey(param_terms_dict, (var, output_idx))
+            param_terms_dict[(var, output_idx)] = zero(T)
         end
-        param_terms_dict[var] += coef
+        param_terms_dict[(var, output_idx)] += coef
     end
 
     return param_terms_dict
 end
 
-# TODO: USED once we update _update_vector_quadratic_constraints!
-# function _delta_parametric_affine_terms(
-#     model,
-#     f::ParametricVectorQuadraticFunction{T},
-# ) where {T}
-#     delta_terms = Dict{Tuple{Int,MOI.VariableIndex},T}()
-
-#     # Handle parameter-variable quadratic terms (px) that become affine (x) when p is updated
-#     for term in f.pv
-#         p_idx_val = p_idx(term.scalar_term.variable_1)
-#         var = term.scalar_term.variable_2
-#         output_idx = term.output_index
-
-#         if haskey(model.updated_parameters, p_idx_val) &&
-#            !isnan(model.updated_parameters[p_idx_val])
-#             old_param_val = model.parameters[p_idx_val]
-#             new_param_val = model.updated_parameters[p_idx_val]
-#             delta_coef =
-#                 term.scalar_term.coefficient * (new_param_val - old_param_val)
-
-#             key = (output_idx, var)
-#             current_delta = get(delta_terms, key, zero(T))
-#             delta_terms[key] = current_delta + delta_coef
-#         end
-#     end
-
-#     # Handle parameter-only affine terms
-#     for term in f.p
-#         p_idx_val = p_idx(term.scalar_term.variable)
-#         output_idx = term.output_index
-
-#         if haskey(model.updated_parameters, p_idx_val) &&
-#            !isnan(model.updated_parameters[p_idx_val])
-#             old_param_val = model.parameters[p_idx_val]
-#             new_param_val = model.updated_parameters[p_idx_val]
-
-#             # This becomes a constant change, not an affine term change
-#             # We'll handle this in the constant update function
-#         end
-#     end
-
-#     return delta_terms
-# end
+function _delta_parametric_constant(
+    model,
+    f::ParametricVectorQuadraticFunction{T},
+) where {T}
+    delta_constants = zeros(T, length(f.current_constant))
+
+    # Handle parameter-only affine terms
+    for term in f.p
+        p_idx_val = p_idx(term.scalar_term.variable)
+        output_idx = term.output_index
+
+        if !isnan(model.updated_parameters[p_idx_val])
+            old_param_val = model.parameters[p_idx_val]
+            new_param_val = model.updated_parameters[p_idx_val]
+            delta_constants[output_idx] +=
+                term.scalar_term.coefficient * (new_param_val - old_param_val)
+        end
+    end
+
+    # Handle parameter-parameter quadratic terms
+    for term in f.pp
+        idx = term.output_index
+        var1 = term.scalar_term.variable_1
+        var2 = term.scalar_term.variable_2
+        p1 = p_idx(var1)
+        p2 = p_idx(var2)
+
+        if !isnan(model.updated_parameters[p1]) ||
+           !isnan(model.updated_parameters[p2])
+            old_val1 = model.parameters[p1]
+            old_val2 = model.parameters[p2]
+            new_val1 =
+                !isnan(model.updated_parameters[p1]) ?
+                model.updated_parameters[p1] : old_val1
+            new_val2 =
+                !isnan(model.updated_parameters[p2]) ?
+                model.updated_parameters[p2] : old_val2
+
+            coef = term.scalar_term.coefficient / (var1 == var2 ? 2 : 1)
+            delta_constants[idx] +=
+                coef * (new_val1 * new_val2 - old_val1 * old_val2)
+        end
+    end
+
+    return delta_constants
+end
+
+function _delta_parametric_affine_terms(
+    model,
+    f::ParametricVectorQuadraticFunction{T},
+) where {T}
+    delta_terms_dict = Dict{Tuple{MOI.VariableIndex,Int},T}()
+    sizehint!(
+        delta_terms_dict,
+        length(vector_quadratic_parameter_variable_terms(f)),
+    )
+
+    # Handle parameter-variable quadratic terms (px) that become affine (x) when p is updated
+    for term in f.pv
+        p_idx_val = p_idx(term.scalar_term.variable_1)
+        var = term.scalar_term.variable_2
+        output_idx = term.output_index
+        if haskey(model.updated_parameters, p_idx_val) &&
+           !isnan(model.updated_parameters[p_idx_val])
+            old_param_val = model.parameters[p_idx_val]
+            new_param_val = model.updated_parameters[p_idx_val]
+            delta_coef =
+                term.scalar_term.coefficient * (new_param_val - old_param_val)
+            base = get(delta_terms_dict, (var, output_idx), zero(T))
+            delta_terms_dict[(var, output_idx)] = base + delta_coef
+        end
+    end
+
+    return delta_terms_dict
+end
 
 function _update_cache!(
     f::ParametricVectorQuadraticFunction{T},
diff --git a/src/update_parameters.jl b/src/update_parameters.jl
index c14db5bd..4f4b08b7 100644
--- a/src/update_parameters.jl
+++ b/src/update_parameters.jl
@@ -160,25 +160,28 @@ function _affine_build_change_and_up_param_func(
     return changes
 end
 
-# TODO
-# This function is currently commented out because there is no vector version
-# of MOI.ScalarCoefficientChange.
-# For the update the entire MOI.ConstraintFunction must be updated
-# function _affine_build_change_and_up_param_func(
-#     pf::ParametricVectorQuadraticFunction{T},
-#     delta_terms,
-# ) where {T}
-#     changes = Vector{MOI.ScalarCoefficientChange}(undef, length(delta_terms))
-#     i = 1
-#     for (var, coef) in delta_terms
-#         base_coef = pf.current_terms_with_p[var]
-#         new_coef = base_coef + coef
-#         pf.current_terms_with_p[var] = new_coef
-#         changes[i] = MOI.ScalarCoefficientChange(var, new_coef)
-#         i += 1
-#     end
-#     return changes
-# end
+function _affine_build_change_and_up_param_func(
+    pf::ParametricVectorQuadraticFunction{T},
+    delta_terms,
+) where {T}
+    for ((var, output_idx), coef) in delta_terms
+        base_coef = pf.current_terms_with_p[(var, output_idx)]
+        new_coef = base_coef + coef
+        pf.current_terms_with_p[(var, output_idx)] = new_coef
+    end
+    new_terms = Dict{MOI.VariableIndex,Vector{Tuple{Int64,T}}}()
+    for ((var, output_idx), coef) in pf.current_terms_with_p
+        if !iszero(coef)
+            base = get!(new_terms, var, Tuple{Int64,T}[])
+            push!(base, (output_idx, coef))
+        end
+    end
+    changes = Vector{MOI.MultirowChange}(undef, length(new_terms))
+    for (i, (var, tuples)) in enumerate(new_terms)
+        changes[i] = MOI.MultirowChange(var, tuples)
+    end
+    return changes
+end
 
 function _update_quadratic_constraints!(
     model::Optimizer,
@@ -290,6 +293,7 @@ function update_parameters!(model::Optimizer)
     _update_affine_constraints!(model)
     _update_vector_affine_constraints!(model)
     _update_quadratic_constraints!(model)
+    _update_vector_quadratic_constraints!(model)
     _update_affine_objective!(model)
     _update_quadratic_objective!(model)
 
@@ -302,17 +306,6 @@ function update_parameters!(model::Optimizer)
         end
     end
 
-    # Different from the above cases, it is not possible to update
-    # vector function in-place because there is no vector version of
-    # MOI.ScalarCoefficientChange.
-    # Therefore, vector quadratic function must be update via
-    # MOI.ConstraintFunction that replaces the entire function.
-    # Consequently, this operation is much slower.
-    # Moreover, to perform the update, it is easier to use the
-    # current_function code, hence, this update must be executed
-    # after the update of the parameters cache (the loop above).
-    _update_vector_quadratic_constraints!(model)
-
     return
 end
 
@@ -330,97 +323,26 @@ function _update_vector_quadratic_constraints!(model::Optimizer)
     return
 end
 
-# TODO: USED once we update _update_vector_quadratic_constraints!
-# function _delta_parametric_constant(
-#     model,
-#     f::ParametricVectorQuadraticFunction{T},
-# ) where {T}
-#     delta_constants = zeros(T, length(f.current_constant))
-
-#     # Handle parameter-only affine terms
-#     for term in f.p
-#         p_idx_val = p_idx(term.scalar_term.variable)
-#         output_idx = term.output_index
-
-#         if !isnan(model.updated_parameters[p_idx_val])
-#             old_param_val = model.parameters[p_idx_val]
-#             new_param_val = model.updated_parameters[p_idx_val]
-#             delta_constants[output_idx] +=
-#                 term.scalar_term.coefficient * (new_param_val - old_param_val)
-#         end
-#     end
-
-#     # Handle parameter-parameter quadratic terms
-#     for term in f.pp
-#         idx = term.output_index
-#         var1 = term.scalar_term.variable_1
-#         var2 = term.scalar_term.variable_2
-#         p1 = p_idx(var1)
-#         p2 = p_idx(var2)
-
-#         if !isnan(model.updated_parameters[p1]) ||
-#            !isnan(model.updated_parameters[p2])
-#             old_val1 = model.parameters[p1]
-#             old_val2 = model.parameters[p2]
-#             new_val1 =
-#                 !isnan(model.updated_parameters[p1]) ?
-#                 model.updated_parameters[p1] : old_val1
-#             new_val2 =
-#                 !isnan(model.updated_parameters[p2]) ?
-#                 model.updated_parameters[p2] : old_val2
-
-#             coef = term.scalar_term.coefficient / (var1 == var2 ? 2 : 1)
-#             delta_constants[idx] +=
-#                 coef * (new_val1 * new_val2 - old_val1 * old_val2)
-#         end
-#     end
-
-#     return delta_constants
-# end
-
-# TODO: Update once MOI.VectorConstantChange is implemented
 function _update_vector_quadratic_constraints!(
-    model::Optimizer,
+    model::Optimizer{T},
     vector_quadratic_constraint_cache_inner::DoubleDictInner{F,S,V},
-) where {F,S,V}
+) where {T,F,S,V}
     for (inner_ci, pf) in vector_quadratic_constraint_cache_inner
-
-        # delta_constants = _delta_parametric_constant(model, pf)
-        # if !iszero(delta_constants)
-        #     pf.current_constant .+= delta_constants
-        #     MOI.modify(
-        #         model.optimizer,
-        #         inner_ci,
-        #         MOI.VectorConstantChange(pf.current_constant),
-        #     )
-        # end
-        # delta_quad_terms = _delta_parametric_affine_terms(model, pf)
-        # if !isempty(delta_quad_terms)
-        #     _quadratic_build_change_and_up_param_func!(pf, delta_quad_terms)
-        #     changes = Vector{MOI.ScalarQuadraticTermChange{T}}()
-        #     for (output_idx, terms) in delta_quad_terms
-        #         for term in terms
-        #             push!(changes, MOI.ScalarQuadraticTermChange(output_idx, term))
-        #         end
-        #     end
-        #     MOI.modify(model.optimizer, inner_ci, changes)
-        # end
-        _update_cache!(pf, model)
-        new_function = _current_function(pf)
-        if _is_vector_affine(new_function)
-            # Build new function if affine
-            new_function = MOI.VectorAffineFunction(
-                new_function.affine_terms,
-                new_function.constants,
+        delta_constants = _delta_parametric_constant(model, pf)
+        if !iszero(delta_constants)
+            pf.current_constant .+= delta_constants
+            MOI.modify(
+                model.optimizer,
+                inner_ci,
+                MOI.VectorConstantChange(pf.current_constant),
             )
         end
-        MOI.set(
-            model.optimizer,
-            MOI.ConstraintFunction(),
-            inner_ci,
-            new_function,
-        )
+        pre_delta_quad_terms = _delta_parametric_affine_terms(model, pf)
+        changes =
+            _affine_build_change_and_up_param_func(pf, pre_delta_quad_terms)
+        for change in changes
+            MOI.modify(model.optimizer, inner_ci, change)
+        end
     end
-
     return
 end
diff --git a/test/jump_tests.jl b/test/jump_tests.jl
index ec3bf6fd..ffaf212a 100644
--- a/test/jump_tests.jl
+++ b/test/jump_tests.jl
@@ -1388,13 +1388,20 @@ function test_jump_psd_cone_with_parameter_pv_v_pv()
         [p * x, (2 * x - 3), p * 3 * x] in
         MOI.PositiveSemidefiniteConeTriangle(2)
     )
+    # which is (p * x) * (p * 3 *x) - (2 * x - 3) ^ 2 >= 0
+    # that simplifies to: p^2 * 3 * x^2 - 4 * x^2 + 12 * x - 9 >= 0
+    # then: (p^2 * 3 - 4) * x^2 + 12 * x - 9 >= 0
+    # for p == 1: -1 * x^2 + 12 * x - 9 >= 0
+    # for p == 3: 23 * x^2 + 12 * x - 9 >= 0
     @objective(model, Min, x)
     @test is_valid(model, con)
     optimize!(model)
     @test value(x) ≈ 0.803845 atol = 1e-5
+    @test value(x) ≈ 6 - 3 * sqrt(3) atol = 1e-5
     set_parameter_value(p, 3.0)
     optimize!(model)
     @test value(x) ≈ 0.416888 atol = 1e-5
+    @test value(x) ≈ (9 * sqrt(3) - 6) / 23 atol = 1e-5
     return delete(model, con)
 end