diff --git a/README.md b/README.md
index 82664b6..01c686c 100644
--- a/README.md
+++ b/README.md
@@ -30,21 +30,21 @@ result = madnlp(model; tol=1e-6)
 #the other input data to generate a structure of named tuples which can then interface with 
 #ExaModels to generate the full model. We do not make any relaxations or decompositions for this problem
 
-model, sc_data, vars, lengths = goc3_model(
+model, cons, vars, lengths, sc_data_array = goc3_model(
     "data/C3E4N00073D1_scenario_303.json", "data/C3E4N00073D1_scenario_303_solution.json"; 
     backend = CUDABackend()
 )
 result = madnlp(model; tol=1e-4)
 
 #Solution from GPU can be used to warm start a CPU solution or vice versa
-model_cpu, sc_data, vars, lengths = goc3_model(
+model, cons, vars, lengths, sc_data_array = goc3_model(
     "data/C3E4N00073D1_scenario_303.json", "data/C3E4N00073D1_scenario_303_solution.json"; 
     result_set = [result, vars]
 )
 result_cpu = ipopt(model_cpu; tol=1e-8)
 
 #Additionally, the SC problem can be evaluated without contingencies
-model, sc_data, vars, lengths = goc3_model(
+model, cons, vars, lengths, sc_data_array = goc3_model(
     "data/C3E4N00073D1_scenario_303.json", "data/C3E4N00073D1_scenario_303_solution.json"; 
     backend = CUDABackend(), include_ctg = false
 )
@@ -99,6 +99,39 @@ result = madnlp(model; tol=1e-6)
 #https://github.com/mit-shin-group/multi-period-opf-data
 ```
 
+### User extension modeling
+ExaModelsPower also supports the user arbitrarily extending any prebuilt models
+```julia
 
+curve = [1, .9, .95]
+# Create vector of NamedTuples elec\_data w/ device data
+untimed_elec_data = [(i = 1, bus = 1, cost = -5000), (i = 2, bus = 2, cost = -2000)]
+Ntime = 3; Nbus = 2
+elec_data = [(;b..., t = t) for b in untimed_elec_data, t in 1:Ntime]
+elec_min = zeros(size(elec_data)); elec_max = fill(50, size(elec_data)); elec_scale = Float64(10)
+
+# User-defined model modifications go here
+function add_electrolyzers(core, vars, cons)
+    # Add new variable to core
+    p_elec = variable(core, size(elec_data, 1), size(elec_data, 2); lvar = elec_min, uvar = elec_max)
+    
+    # Objectives are additive. Add secondary objective
+    o2 = objective(core, e.cost*p_elec[e.i, e.t] for e in elec_data)
+    
+    # Add electrolyzer load to power balance for each bus
+    c_elec_power_balance = constraint!(core, cons.c_active_power_balance, e.bus + Nbus*(e.t-1) => p_elec[e.i, e.t] for e in elec_data)
+    
+    # Ramping limit over time
+    c_elec_ramp = constraint(core, p_elec[e.i, e.t] - p_elec[e.i, e.t - 1] for e in elec_data[:, 2:Ntime]; lcon = fill(-elec_scale, size(elec_data[:, 2:Ntime])), ucon = fill(elec_scale, size(elec_data[:, 2:Ntime])))
+    
+    # Set initial electrolyzer power to 0
+    c_elec_ramp_init = constraint(core, p_elec[e.i, e.t] for e in elec_data[:, 1];)
+
+    # Return new variables and constraints to be tracked
+    return ((p_elec=p_elec,), (c_elec_ramp=c_elec_ramp, c_elec_ramp_init=c_elec_ramp_init))
+end
+# Generate model
+model, vars, cons = mpopf_model("pglib_opf_case3_lmbd.m", curve; user_callback = add_electrolyzers) # user_callback function added after initial mpopf model is constructed
+```
 
 
diff --git a/src/mpopf.jl b/src/mpopf.jl
index b9289ec..28fa4bd 100644
--- a/src/mpopf.jl
+++ b/src/mpopf.jl
@@ -234,8 +234,11 @@ function build_mpopf(data, Nbus, N, form, user_callback; backend = nothing, T =
         vars, cons = add_piecewise_cons(core, data, N, vars, cons, storage_complementarity_constraint)
     end
 
-    vars, cons = user_callback(core, vars, cons)
-    model = ExaModel(core; kwargs...)
+    vars2, cons2 = user_callback(core, vars, cons)
+    model =ExaModel(core; kwargs...)
+
+    vars = (;vars..., vars2...)
+    cons = (;cons..., cons2...)
     return model, vars, cons
 end
 
@@ -251,8 +254,11 @@ function build_mpopf(data, Nbus, N, discharge_func::Function, form, user_callbac
         vars, cons = add_smooth_cons(core, data, N, vars, cons, discharge_func)
     end
 
-    vars, cons = user_callback(core, vars, cons)
-    model = ExaModel(core; kwargs...)
+    vars2, cons2 = user_callback(core, vars, cons)
+    model =ExaModel(core; kwargs...)
+
+    vars = (;vars..., vars2...)
+    cons = (;cons..., cons2...)
     return model, vars, cons
 end
 
diff --git a/src/opf.jl b/src/opf.jl
index 0a298cd..3932f54 100644
--- a/src/opf.jl
+++ b/src/opf.jl
@@ -1,5 +1,5 @@
 function dummy_extension(core, vars, cons)
-    return vars, cons
+    return (;), (;)
 end
 
 function build_polar_opf(data, user_callback; backend = nothing, T=Float64, kwargs...)
@@ -87,9 +87,12 @@ function build_polar_opf(data, user_callback; backend = nothing, T=Float64, kwar
         c_to_thermal_limit = c_to_thermal_limit
     )
 
-    vars, cons = user_callback(core, vars, cons)
+    vars2, cons2 = user_callback(core, vars, cons)
     model =ExaModel(core; kwargs...)
 
+    vars = (;vars..., vars2...)
+    cons = (;cons..., cons2...)
+
     return model, vars, cons
 end
 
@@ -182,8 +185,11 @@ function build_rect_opf(data, user_callback; backend = nothing, T=Float64, kwarg
         c_voltage_magnitude = c_voltage_magnitude
     )
 
-    vars, cons = user_callback(core, vars, cons)
-    model = ExaModel(core; kwargs...)
+    vars2, cons2 = user_callback(core, vars, cons)
+    model =ExaModel(core; kwargs...)
+
+    vars = (;vars..., vars2...)
+    cons = (;cons..., cons2...)
 
     return model, vars, cons
 end
diff --git a/src/scopf.jl b/src/scopf.jl
index 2855793..d63b799 100644
--- a/src/scopf.jl
+++ b/src/scopf.jl
@@ -4,6 +4,7 @@ function goc3_model(
     T = Float64,
     include_ctg = true, #contingencies will be specified in input data
     result_set = [],
+    user_callback = dummy_extension,
     kwargs...
         )
 
@@ -241,10 +242,11 @@ function goc3_model(
     end
 
     #4.2.1 Bus power mismatch and penalized mismatch definitions
-    c_11 = constraint(core, p_it_plus[b.i, b.t] - p_it[b.i, b.t] for b in sc_data.busarray; ucon = fill(Inf, size(sc_data.busarray)))
-    c_12 = constraint(core, p_it_plus[b.i, b.t] + p_it[b.i, b.t] for b in sc_data.busarray; ucon = fill(Inf, size(sc_data.busarray)))
-    c_13 = constraint(core, q_it_plus[b.i, b.t] - q_it[b.i, b.t] for b in sc_data.busarray; ucon = fill(Inf, size(sc_data.busarray)))
-    c_14 = constraint(core, q_it_plus[b.i, b.t] + q_it[b.i, b.t] for b in sc_data.busarray; ucon = fill(Inf, size(sc_data.busarray)))
+    c11 = constraint(core, p_it_plus[b.i, b.t] - p_it[b.i, b.t] for b in sc_data.busarray; ucon = fill(Inf, size(sc_data.busarray)))
+    c12 = constraint(core, p_it_plus[b.i, b.t] + p_it[b.i, b.t] for b in sc_data.busarray; ucon = fill(Inf, size(sc_data.busarray)))
+    c13 = constraint(core, q_it_plus[b.i, b.t] - q_it[b.i, b.t] for b in sc_data.busarray; ucon = fill(Inf, size(sc_data.busarray)))
+    c14 = constraint(core, q_it_plus[b.i, b.t] + q_it[b.i, b.t] for b in sc_data.busarray; ucon = fill(Inf, size(sc_data.busarray)))
+    cons = (c11=c11, c12=c12, c13=c13, c14=c14)
     #4.2.2 Bus pwoer mismatch penalty
     c15 = constraint(core, z_it_p[b.i, b.t] - b.dt*sum(sc_data.c_p)*p_it_plus[b.i, b.t] for b in sc_data.busarray)
     c16 = constraint(core, z_it_q[b.i, b.t] - b.dt*sum(sc_data.c_q)*q_it_plus[b.i, b.t] for b in sc_data.busarray)
@@ -254,6 +256,7 @@ function goc3_model(
     c17_cs = constraint!(core, c17, cs.bus + I*(cs.t-1) => -p_jt_cs[cs.j_cs, cs.t] for cs in sc_data.csarray)
     c17_sh = constraint!(core, c17, sh.bus + I*(sh.t-1) => -p_jt_sh[sh.j_sh, sh.t] for sh in sc_data.shuntarray)
     
+    
     #Reminder: fr and to split for ln, xf, and dc
     c17_fr_ln = constraint!(core, c17, ln.fr_bus + I*(ln.t-1) => -p_jt_fr_ln[ln.j_ln, ln.t] for ln in sc_data.aclbrancharray)
     c17_fr_xf = constraint!(core, c17, xf.fr_bus + I*(xf.t-1) => -p_jt_fr_xf[xf.j_xf, xf.t] for xf in sc_data.acxbrancharray)
@@ -272,7 +275,6 @@ function goc3_model(
     c18_to_ln = constraint!(core, c18, ln.to_bus + I*(ln.t-1) => -q_jt_to_ln[ln.j_ln, ln.t] for ln in sc_data.aclbrancharray)
     c18_to_xf = constraint!(core, c18, xf.to_bus + I*(xf.t-1) => -q_jt_to_xf[xf.j_xf, xf.t] for xf in sc_data.acxbrancharray)
     c18_to_dc = constraint!(core, c18, dc.to_bus + I*(dc.t-1) => -q_jt_to_dc[dc.j_dc, dc.t] for dc in sc_data.dclinearray)
-
     
     #4.3.2 Reserve shortfall penalties
     c28 = constraint(core, z_nt_rgu[n.n_p, n.t] - n.dt*n.c_rgu*p_nt_rgu_plus[n.n_p, n.t] for n in sc_data.preservearray)
@@ -283,7 +285,6 @@ function goc3_model(
     c33 = constraint(core, z_nt_rrd[n.n_p, n.t] - n.dt*n.c_rrd*p_nt_rrd_plus[n.n_p, n.t] for n in sc_data.preservearray)
     c34 = constraint(core, z_nt_qru[n.n_q, n.t] - n.dt*n.c_qru*q_nt_qru_plus[n.n_q, n.t] for n in sc_data.qreservearray)
     c35 = constraint(core, z_nt_qrd[n.n_q, n.t] - n.dt*n.c_qrd*q_nt_qrd_plus[n.n_q, n.t] for n in sc_data.qreservearray)
-
     
     #4.3.3 Reserve requirements
     c36 = constraint(core, p_nt_rgu_req[n.n_p, n.t]/n.σ_rgu for n in sc_data.preservearray)
@@ -477,7 +478,6 @@ function goc3_model(
 
     #4.8.3 AC branch flows
     
-    
     c148_ln = constraint(core, -p_jt_fr_ln[ln.j_ln, ln.t] + ln.u_on*((ln.g_sr+ln.g_fr)*v_it[ln.fr_bus, ln.t]^2 + (-ln.g_sr*cos(θ_it[ln.fr_bus, ln.t] - θ_it[ln.to_bus, ln.t])
      - ln.b_sr*sin(θ_it[ln.fr_bus, ln.t] - θ_it[ln.to_bus, ln.t]))*v_it[ln.fr_bus, ln.t]*v_it[ln.to_bus, ln.t]) for ln in sc_data.aclbrancharray)
     c148_xf = constraint(core, -p_jt_fr_xf[xf.j_xf, xf.t] + xf.u_on*((xf.g_sr+xf.g_fr)*v_it[xf.fr_bus, xf.t]^2/(τ_jt_xf[xf.j_xf, xf.t]^2) + (-xf.g_sr*cos(θ_it[xf.fr_bus, xf.t] - θ_it[xf.to_bus, xf.t] - φ_jt_xf[xf.j_xf, xf.t])
@@ -494,7 +494,6 @@ function goc3_model(
     + ln.g_sr*sin(θ_it[ln.fr_bus, ln.t] - θ_it[ln.to_bus, ln.t]))*v_it[ln.fr_bus, ln.t]*v_it[ln.to_bus, ln.t]) for ln in sc_data.aclbrancharray)
     c151_xf = constraint(core, -q_jt_to_xf[xf.j_xf, xf.t] + xf.u_on*((-xf.b_sr-xf.b_to-xf.b_ch/2)*v_it[xf.to_bus, xf.t]^2 + (xf.b_sr*cos(θ_it[xf.fr_bus, xf.t] - θ_it[xf.to_bus, xf.t] - φ_jt_xf[xf.j_xf, xf.t]) 
     + xf.g_sr*sin(θ_it[xf.fr_bus, xf.t] - θ_it[xf.to_bus, xf.t] - φ_jt_xf[xf.j_xf, xf.t]))*v_it[xf.fr_bus, xf.t]*v_it[xf.to_bus, xf.t]/τ_jt_xf[xf.j_xf, xf.t]) for xf in sc_data.acxbrancharray)
-    
 
     #4.8.4 DC lines
     c156 = constraint(core, p_jt_fr_dc[dc.j_dc, dc.t] + p_jt_to_dc[dc.j_dc, dc.t] for dc in sc_data.dclinearray)
@@ -535,11 +534,10 @@ function goc3_model(
         c163_sh = constraint!(core, c163, sh.bus + I*(sh.t-1) + I*L_T*(sh.k-1) => -p_jt_sh[sh.j_sh, sh.t] for sh in sc_data.k_shuntarray)
         c163_dc_fr = constraint!(core, c163, dc.fr_bus + I*(dc.t-1) + I*L_T*(dc.ctg-1) => -p_jt_fr_dc[dc.j_dc, dc.t] for dc in sc_data.jtk_dc_flattened)
         c163_dc_to = constraint!(core, c163, dc.to_bus + I*(dc.t-1) + I*L_T*(dc.ctg-1) => -p_jt_fr_dc[dc.j_dc, dc.t] for dc in sc_data.jtk_dc_flattened)
+        
     end
 
 
-    model = ExaModel(core; kwargs...)
-
     vars = (b_jt_sh = b_jt_sh,
             g_jt_sh = g_jt_sh,
             #Split e_w_plus into separate sets for W_en_min and W_en_max ad for pr, cs
@@ -657,6 +655,128 @@ function goc3_model(
             τ_jt_xf = τ_jt_xf,
             φ_jt_xf = φ_jt_xf)
 
+        cons = (
+            c11 = c11,
+            c12 = c12,
+            c13 = c13,
+            c14 = c14,
+            c15 = c15,
+            c16 = c16,
+            c17 = c17,
+            c18 = c18,
+            c28 = c28,
+            c29 = c29,
+            c30 = c30,
+            c31 = c31,
+            c32 = c32,
+            c33 = c33,
+            c34 = c34,
+            c35 = c35,
+            c36 = c36,
+            c37 = c37,
+            cmax38 = cmax38,
+            c38 = c38,
+            c39 = c39,
+            c40 = c40,
+            c41 = c41,
+            c42 = c42,
+            c43 = c43,
+            c44 = c44,
+            c45 = c45,
+            c46 = c46,
+            c47 = c47,
+            c68_pr = c68_pr,
+            c68_cs = c68_cs,
+            c69_pr = c69_pr,
+            c69_cs = c69_cs,
+            c70_pr = c70_pr,
+            c70_cs = c70_cs,
+            c71_pr = c71_pr,
+            c71_cs = c71_cs,
+            c72_pr = c72_pr,
+            c72_cs = c72_cs,
+            c73_pr = c73_pr,
+            c73_cs = c73_cs,
+            c74_pr = c74_pr,
+            c74_cs = c74_cs,
+            c75_pr = c75_pr,
+            c75_cs = c75_cs,
+            c76_pr = c76_pr,
+            c76_cs = c76_cs,
+            c78_pr = c78_pr,
+            c78_cs = c78_cs,
+            c79_pr = c79_pr,
+            c79_cs = c79_cs,
+            c94_pr = c94_pr,
+            c94_cs = c94_cs,
+            c95_pr = c95_pr,
+            c95_cs = c95_cs,
+            c96_pr = c96_pr,
+            c96_cs = c96_cs,
+            c97_pr = c97_pr,
+            c97_cs = c97_cs,
+            c98_pr = c98_pr,
+            c98_cs = c98_cs,
+            c99_pr = c99_pr,
+            c99_cs = c99_cs,
+            c100_pr = c100_pr,
+            c100_cs = c100_cs,
+            c101_pr = c101_pr,
+            c102_pr = c102_pr,
+            c102_cs = c102_cs,
+            c103_pr = c103_pr,
+            c104_pr = c104_pr,
+            c104_cs = c104_cs,
+            c105_cs = c105_cs,
+            c109 = c109,
+            c110 = c110,
+            c111 = c111,
+            c112 = c112,
+            c113 = c113,
+            c114 = c114,
+            c115 = c115,
+            c116 = c116,
+            c117 = c117,
+            c118 = c118,
+            c119 = c119,
+            c120 = c120,
+            c121 = c121,
+            c122 = c122,
+            c123 = c123,
+            c124 = c124,
+            c125 = c125,
+            c126 = c126,
+            c127 = c127,
+            c128 = c128,
+            c130_pr = c130_pr,
+            c130_cs = c130_cs,
+            c131_pr = c131_pr,
+            c131_cs = c131_cs,
+            c132 = c132,
+            c133 = c133,
+            c134 = c134,
+            c135 = c135,
+            c139_ln = c139_ln,
+            c139_xf = c139_xf,
+            c140_ln = c140_ln,
+            c140_xf = c140_xf,
+            c141_ln = c141_ln,
+            c141_xf = c141_xf,
+            c144 = c144,
+            c145 = c145,
+            c146 = c146,
+            c147 = c147,
+            c148_ln = c148_ln,
+            c148_xf = c148_xf,
+            c149_ln = c149_ln,
+            c149_xf = c149_xf,
+            c150_ln = c150_ln,
+            c150_xf = c150_xf,
+            c151_ln = c151_ln,
+            c151_xf = c151_xf,
+            c156 = c156
+        )
+
         if include_ctg
             vars = (;vars...,
             s_jtk_plus_ln = s_jtk_plus_ln,
@@ -670,11 +790,30 @@ function goc3_model(
             z_tk_ctg = z_tk_ctg,
             z_t_ctg_min = z_t_ctg_min,
             z_t_ctg_avg = z_t_ctg_avg)
+            cons = (;cons...,
+                    c4=c4, 
+                    c5=c5, 
+                    c10=c10, 
+                    c158_ln=c158_ln, 
+                    c158_xf=c158_xf, 
+                    c159_ln=c159_ln, 
+                    c159_xf=c159_xf, 
+                    c160_ln=c160_ln, 
+                    c160_xf=c160_xf, 
+                    c161_ln=c161_ln, 
+                    c161_xf=c161_xf, 
+                    c162=c162, 
+                    c163=c163)
         end
 
+        vars2, cons2 = user_callback(core, vars, cons)
+        model =ExaModel(core; kwargs...)
+
+        vars = (;vars..., vars2...)
+        cons = (;cons..., cons2...)
    
     @info "built model"
-    return model, sc_data_array, vars, lengths
+    return model, cons, vars, lengths, sc_data_array
 
 end
 
diff --git a/test/opf_tests.jl b/test/opf_tests.jl
index d3aaf72..9d30439 100644
--- a/test/opf_tests.jl
+++ b/test/opf_tests.jl
@@ -69,6 +69,6 @@ end
 function sc_tests(filename, backend, T)
     uc_filename = filename*"_solution.json"
     filename = filename*".json"
-    model, sc_data_array, vars, lengths = ExaModelsPower.goc3_model(filename, uc_filename; backend=backend, T=T)
+    model, cons, vars, lengths, sc_data_array = ExaModelsPower.goc3_model(filename, uc_filename; backend=backend, T=T)
     result = madnlp(model; max_iter=5, tol=1e-2)
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 997a6e7..76f07ad 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -83,8 +83,8 @@ function add_electrolyzers(core, vars, cons)
     length(elec_data)), -elec_scale),
     ucon = fill!(similar(elec_data, Float64,
     length(elec_data)), elec_scale))
-    vars = (;vars..., p_elec=p_elec)
-    cons = (;cons..., c_elec_ramp=c_elec_ramp)
+    vars = (p_elec=p_elec,)
+    cons = (c_elec_ramp=c_elec_ramp,)
     return vars, cons
 end
 
@@ -224,7 +224,7 @@ function runtests()
 
             @testset "User callback, $(T), $(backend)" begin
                 model, vars, cons = mpopf_model(
-                    "../data/pglib_opf_case3_lmbd.m", elec_curve;
+                    "../data/pglib_opf_case3_lmbd_mod.m", elec_curve;
                     user_callback = add_electrolyzers, T=T, backend=backend)
             end