diff --git a/src/model/build.jl b/src/model/build.jl
index a9bf454..fa0ec0c 100644
--- a/src/model/build.jl
+++ b/src/model/build.jl
@@ -16,7 +16,6 @@ function build_model(instance::Instance; optimizer, variable_names::Bool = false
     E_out = Dict(src => [] for src in plants ∪ centers)
 
     function push_edge!(src, dst, m)
-        @show src.name, dst.name, m.name
         push!(E, (src, dst, m))
         push!(E_out[src], (dst, m))
         push!(E_in[dst], (src, m))
@@ -35,7 +34,6 @@ function build_model(instance::Instance; optimizer, variable_names::Bool = false
 
             # Plant to center
             for c in centers
-                @show m.name, p1.name, c.name, m == c.input
                 m == c.input || continue
                 push_edge!(p1, c, m)
             end
@@ -63,7 +61,6 @@ function build_model(instance::Instance; optimizer, variable_names::Bool = false
     for (p1, p2, m) in E
         d = _calculate_distance(p1.latitude, p1.longitude, p2.latitude, p2.longitude)
         distances[p1, p2, m] = d
-        @show p1.name, p2.name, m.name, d
     end
 
     # Decision variables
@@ -102,11 +99,20 @@ function build_model(instance::Instance; optimizer, variable_names::Bool = false
         z_disp[c.name, m.name, t] = @variable(model, lower_bound = 0)
     end
 
-    # Total plant input
+    # Total plant/center input
     z_input = _init(model, :z_input)
     for p in plants, t in T
         z_input[p.name, t] = @variable(model, lower_bound = 0)
     end
+    for c in centers, t in T
+        z_input[c.name, t] = @variable(model, lower_bound = 0)
+    end
+
+    # Total amount collected by the center
+    z_collected = _init(model, :z_collected)
+    for c in centers, m in c.outputs, t in T
+        z_collected[c.name, m.name, t] = @variable(model, lower_bound = 0)
+    end
 
 
     # Objective function
@@ -115,47 +121,55 @@ function build_model(instance::Instance; optimizer, variable_names::Bool = false
 
     # Transportation cost
     for (p1, p2, m) in E, t in T
-        obj += distances[p1, p2, m] * y[p1.name, p2.name, m.name, t]
+        add_to_expression!(obj, distances[p1, p2, m], y[p1.name, p2.name, m.name, t])
     end
 
     # Center: Revenue
     for c in centers, (p, m) in E_in[c], t in T
-        obj += c.revenue[t] * y[p.name, c.name, m.name, t]
+        add_to_expression!(obj, c.revenue[t], y[p.name, c.name, m.name, t])
     end
 
     # Center: Collection cost
     for c in centers, (p, m) in E_out[c], t in T
-        obj += c.collection_cost[m][t] * y[c.name, p.name, m.name, t]
+        add_to_expression!(obj, c.collection_cost[m][t], y[c.name, p.name, m.name, t])
     end
 
     # Center: Disposal cost
     for c in centers, m in c.outputs, t in T
-        obj += c.disposal_cost[m][t] * z_disp[c.name, m.name, t]
+        add_to_expression!(obj, c.disposal_cost[m][t], z_disp[c.name, m.name, t])
     end
 
     # Center: Operating cost
     for c in centers, t in T
-        obj += c.operating_cost[t]
+        add_to_expression!(obj, c.operating_cost[t])
     end
 
     # Plants: Disposal cost
     for p in plants, m in keys(p.output), t in T
-        obj += p.disposal_cost[m][t] * z_disp[p.name, m.name, t]
+        add_to_expression!(obj, p.disposal_cost[m][t], z_disp[p.name, m.name, t])
     end
 
     # Plants: Opening cost
     for p in plants, t in T
-        obj += p.capacities[1].opening_cost[t] * (x[p.name, t] - x[p.name, t-1])
+        add_to_expression!(
+            obj,
+            p.capacities[1].opening_cost[t],
+            (x[p.name, t] - x[p.name, t-1]),
+        )
     end
 
     # Plants: Fixed operating cost
     for p in plants, t in T
-        obj += p.capacities[1].fix_operating_cost[t] * x[p.name, t]
+        add_to_expression!(obj, p.capacities[1].fix_operating_cost[t], x[p.name, t])
     end
 
     # Plants: Variable operating cost
     for p in plants, (src, m) in E_in[p], t in T
-        obj += p.capacities[1].var_operating_cost[t] * y[src.name, p.name, m.name, t]
+        add_to_expression!(
+            obj,
+            p.capacities[1].var_operating_cost[t],
+            y[src.name, p.name, m.name, t],
+        )
     end
 
     @objective(model, Min, obj)
@@ -206,31 +220,75 @@ function build_model(instance::Instance; optimizer, variable_names::Bool = false
     # Plants: Capacity limit
     eq_capacity = _init(model, :eq_capacity)
     for p in plants, t in T
-        eq_capacity[p.name, t] = @constraint(
-            model,
-            z_input[p.name, t] <= p.capacities[1].size * x[p.name, t]
-        )
+        eq_capacity[p.name, t] =
+            @constraint(model, z_input[p.name, t] <= p.capacities[1].size * x[p.name, t])
     end
 
     # Plants: Disposal limit
     eq_disposal_limit = _init(model, :eq_disposal_limit)
     for p in plants, m in keys(p.output), t in T
         isfinite(p.disposal_limit[m][t]) || continue
-        eq_disposal_limit[p.name, m.name, t] = @constraint(
-            model,
-            z_disp[p.name, m.name, t] <= p.disposal_limit[m][t]
-        )
+        eq_disposal_limit[p.name, m.name, t] =
+            @constraint(model, z_disp[p.name, m.name, t] <= p.disposal_limit[m][t])
     end
 
     # Plants: Plant remains open
     eq_keep_open = _init(model, :eq_keep_open)
     for p in plants, t in T
-        eq_keep_open[p.name, t] = @constraint(
+        eq_keep_open[p.name, t] = @constraint(model, x[p.name, t] >= x[p.name, t-1])
+    end
+
+    # Plants: Building period
+    eq_building_period = _init(model, :eq_building_period)
+    for p in plants, t in T
+        if t ∉ instance.building_period
+            eq_building_period[p.name, t] = @constraint(model, x[p.name, t] == 0)
+        end
+    end
+
+    # Centers: Definition of total center input
+    eq_z_input = _init(model, :eq_z_input)
+    for c in centers, t in T
+        eq_z_input[c.name, t] = @constraint(
             model,
-            x[p.name, t] >= x[p.name, t-1]
+            z_input[c.name, t] ==
+            sum(y[src.name, c.name, m.name, t] for (src, m) in E_in[c])
         )
     end
 
+    # Centers: Calculate amount collected
+    eq_z_collected = _init(model, :eq_z_collected)
+    for c in centers, m in c.outputs, t in T
+        M = length(c.var_output[m])
+        eq_z_collected[c.name, m.name, t] = @constraint(
+            model,
+            z_collected[c.name, m.name, t] ==
+            sum(
+                z_input[c.name, t-offset] * c.var_output[m][offset+1] for
+                offset = 0:min(M - 1, t - 1)
+            ) + c.fixed_output[m][t]
+        )
+    end
+
+    # Centers: Collected products must be disposed or sent
+    eq_balance = _init(model, :eq_balance)
+    for c in centers, m in c.outputs, t in T
+        eq_balance[c.name, m.name, t] = @constraint(
+            model,
+            z_collected[c.name, m.name, t] ==
+            sum(y[c.name, dst.name, m.name, t] for (dst, m2) in E_out[c] if m == m2) +
+            z_disp[c.name, m.name, t]
+        )
+    end
+
+    # Centers: Disposal limit
+    eq_disposal_limit = _init(model, :eq_disposal_limit)
+    for c in centers, m in c.outputs, t in T
+        isfinite(c.disposal_limit[m][t]) || continue
+        eq_disposal_limit[c.name, m.name, t] =
+            @constraint(model, z_disp[c.name, m.name, t] <= c.disposal_limit[m][t])
+    end
+
     if variable_names
         _set_names!(model)
     end
diff --git a/test/fixtures/simple.json b/test/fixtures/simple.json
index bb104c5..c5b7f1e 100644
--- a/test/fixtures/simple.json
+++ b/test/fixtures/simple.json
@@ -49,8 +49,8 @@
         "P3": [20, 10, 0, 0]
       },
       "variable output (tonne/tonne)": {
-        "P2": [0.12, 0.25, 0.12, 0.0],
-        "P3": [0.25, 0.25, 0.25, 0.0]
+        "P2": [0.20, 0.25, 0.12],
+        "P3": [0.25, 0.25, 0.25]
       },
       "revenue ($/tonne)": 12.0,
       "collection cost ($/tonne)": {
diff --git a/test/src/RELOGT.jl b/test/src/RELOGT.jl
index a90e971..53a92ad 100644
--- a/test/src/RELOGT.jl
+++ b/test/src/RELOGT.jl
@@ -18,7 +18,8 @@ function runtests()
     @testset "RELOG" begin
         instance_parse_test_1()
         instance_parse_test_2()
-        model_build_test()
+        model_build_test_1()
+        model_build_test_2()
         model_dist_test()
     end
 end
diff --git a/test/src/instance/parse_test.jl b/test/src/instance/parse_test.jl
index 270556e..07697c4 100644
--- a/test/src/instance/parse_test.jl
+++ b/test/src/instance/parse_test.jl
@@ -31,8 +31,7 @@ function instance_parse_test_1()
     @test c1.input === p1
     @test c1.outputs == [p2, p3]
     @test c1.fixed_output == Dict(p2 => [100, 50, 0, 0], p3 => [20, 10, 0, 0])
-    @test c1.var_output ==
-          Dict(p2 => [0.12, 0.25, 0.12, 0.0], p3 => [0.25, 0.25, 0.25, 0.0])
+    @test c1.var_output == Dict(p2 => [0.2, 0.25, 0.12], p3 => [0.25, 0.25, 0.25])
     @test c1.revenue == [12.0, 12.0, 12.0, 12.0]
     @test c1.operating_cost == [150.0, 150.0, 150.0, 150.0]
     @test c1.disposal_limit == Dict(p2 => [0, 0, 0, 0], p3 => [Inf, Inf, Inf, Inf])
diff --git a/test/src/model/build_test.jl b/test/src/model/build_test.jl
index a630fce..85764e0 100644
--- a/test/src/model/build_test.jl
+++ b/test/src/model/build_test.jl
@@ -3,7 +3,7 @@ using Test
 using HiGHS
 using JuMP
 
-function model_build_test()
+function model_build_test_1()
     instance = RELOG.parsefile(fixture("simple.json"))
     model = RELOG.build_model(instance, optimizer = HiGHS.Optimizer, variable_names = true)
     y = model[:y]
@@ -50,34 +50,70 @@ function model_build_test()
 
     # Plants: Must meet input mix
     @test repr(model[:eq_input_mix]["L1", "P1", 1]) ==
-        "eq_input_mix[L1,P1,1] : y[C2,L1,P1,1] - 0.953 z_input[L1,1] = 0"
+          "eq_input_mix[L1,P1,1] : y[C2,L1,P1,1] - 0.953 z_input[L1,1] = 0"
     @test repr(model[:eq_input_mix]["L1", "P2", 1]) ==
-        "eq_input_mix[L1,P2,1] : y[C1,L1,P2,1] - 0.047 z_input[L1,1] = 0"
+          "eq_input_mix[L1,P2,1] : y[C1,L1,P2,1] - 0.047 z_input[L1,1] = 0"
 
     # Plants: Calculate amount produced
     @test repr(model[:eq_z_prod]["L1", "P3", 1]) ==
-        "eq_z_prod[L1,P3,1] : z_prod[L1,P3,1] - 0.25 z_input[L1,1] = 0"
+          "eq_z_prod[L1,P3,1] : z_prod[L1,P3,1] - 0.25 z_input[L1,1] = 0"
     @test repr(model[:eq_z_prod]["L1", "P4", 1]) ==
-        "eq_z_prod[L1,P4,1] : z_prod[L1,P4,1] - 0.12 z_input[L1,1] = 0"
+          "eq_z_prod[L1,P4,1] : z_prod[L1,P4,1] - 0.12 z_input[L1,1] = 0"
 
     # Plants: Produced material must be sent or disposed
     @test repr(model[:eq_balance]["L1", "P3", 1]) ==
-        "eq_balance[L1,P3,1] : z_prod[L1,P3,1] - z_disp[L1,P3,1] = 0"
+          "eq_balance[L1,P3,1] : z_prod[L1,P3,1] - z_disp[L1,P3,1] = 0"
     @test repr(model[:eq_balance]["L1", "P4", 1]) ==
-        "eq_balance[L1,P4,1] : -y[L1,C3,P4,1] + z_prod[L1,P4,1] - z_disp[L1,P4,1] = 0"
+          "eq_balance[L1,P4,1] : -y[L1,C3,P4,1] + z_prod[L1,P4,1] - z_disp[L1,P4,1] = 0"
 
     # Plants: Capacity limit
     @test repr(model[:eq_capacity]["L1", 1]) ==
-        "eq_capacity[L1,1] : -100 x[L1,1] + z_input[L1,1] ≤ 0"
+          "eq_capacity[L1,1] : -100 x[L1,1] + z_input[L1,1] ≤ 0"
 
     # Plants: Disposal limit
     @test repr(model[:eq_disposal_limit]["L1", "P4", 1]) ==
-        "eq_disposal_limit[L1,P4,1] : z_disp[L1,P4,1] ≤ 1000"
+          "eq_disposal_limit[L1,P4,1] : z_disp[L1,P4,1] ≤ 1000"
     @test ("L1", "P3", 1) ∉ keys(model[:eq_disposal_limit])
 
     # Plants: Plant remains open
     @test repr(model[:eq_keep_open]["L1", 4]) ==
-        "eq_keep_open[L1,4] : -x[L1,3] + x[L1,4] ≥ 0"
+          "eq_keep_open[L1,4] : -x[L1,3] + x[L1,4] ≥ 0"
     @test repr(model[:eq_keep_open]["L1", 1]) == "eq_keep_open[L1,1] : x[L1,1] ≥ 0"
 
+    # Plants: Building period
+    @test ("L1", 1) ∉ keys(model[:eq_building_period])
+    @test repr(model[:eq_building_period]["L1", 2]) ==
+          "eq_building_period[L1,2] : x[L1,2] = 0"
+
+    # Centers: Definition of total center input
+    @test repr(model[:eq_z_input]["C1", 1]) ==
+          "eq_z_input[C1,1] : -y[C2,C1,P1,1] + z_input[C1,1] = 0"
+
+    # Centers: Calculate amount collected
+    @test repr(model[:eq_z_collected]["C1", "P2", 1]) ==
+          "eq_z_collected[C1,P2,1] : -0.2 z_input[C1,1] + z_collected[C1,P2,1] = 100"
+    @test repr(model[:eq_z_collected]["C1", "P2", 2]) ==
+          "eq_z_collected[C1,P2,2] : -0.25 z_input[C1,1] - 0.2 z_input[C1,2] + z_collected[C1,P2,2] = 50"
+    @test repr(model[:eq_z_collected]["C1", "P2", 3]) ==
+          "eq_z_collected[C1,P2,3] : -0.12 z_input[C1,1] - 0.25 z_input[C1,2] - 0.2 z_input[C1,3] + z_collected[C1,P2,3] = 0"
+    @test repr(model[:eq_z_collected]["C1", "P2", 4]) ==
+          "eq_z_collected[C1,P2,4] : -0.12 z_input[C1,2] - 0.25 z_input[C1,3] - 0.2 z_input[C1,4] + z_collected[C1,P2,4] = 0"
+
+    # Centers: Collected products must be disposed or sent
+    @test repr(model[:eq_balance]["C1", "P2", 1]) ==
+          "eq_balance[C1,P2,1] : -y[C1,L1,P2,1] - z_disp[C1,P2,1] + z_collected[C1,P2,1] = 0"
+    @test repr(model[:eq_balance]["C1", "P3", 1]) ==
+          "eq_balance[C1,P3,1] : -z_disp[C1,P3,1] + z_collected[C1,P3,1] = 0"
+
+    # Centers: Disposal limit
+    @test repr(model[:eq_disposal_limit]["C1", "P2", 1]) ==
+          "eq_disposal_limit[C1,P2,1] : z_disp[C1,P2,1] ≤ 0"
+    @test ("C1", "P3", 1) ∉ keys(model[:eq_disposal_limit])
+end
+
+
+function model_build_test_2()
+    instance = RELOG.parsefile(fixture("boat_example.json"))
+    model = RELOG.build_model(instance, optimizer = HiGHS.Optimizer)
+    optimize!(model) 
 end