diff --git a/REQUIRE b/REQUIRE
index 137767a..57f9d3b 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -1 +1,6 @@
 julia 0.6
+DataFrames
+Iterators
+TransitNetworks
+JuMP
+Gurobi
\ No newline at end of file
diff --git a/src/NetworkDesignModels.jl b/src/NetworkDesignModels.jl
index 5c8d7bd..08b6d45 100644
--- a/src/NetworkDesignModels.jl
+++ b/src/NetworkDesignModels.jl
@@ -1,5 +1,13 @@
 module NetworkDesignModels
+    
+    export DesignModel
 
-# package code goes here
+    import DataFrames, Iterators, JuMP, Gurobi
+    using TransitNetworks
 
-end # module
+    include("model/designmodel.jl")
+    include("model/variables.jl")
+    include("model/constraints.jl")
+    include("model/objective.jl")
+
+end
diff --git a/src/model/constraints.jl b/src/model/constraints.jl
new file mode 100644
index 0000000..3a5ed90
--- /dev/null
+++ b/src/model/constraints.jl
@@ -0,0 +1,49 @@
+function routeopenconstraint(dm::DesignModel)
+    validod = vec([isassigned(dm.np.stage, u,v) 
+        for u in 1:dm.np.nstations, v in 1:dm.np.nstations])
+    stages = unique(vcat(dm.np.stage[validod]...))
+    JuMP.@constraint(dm.model, 
+        [r=1:length(stages), s=stages[r]], 
+        dm.open[stages[r]] <= dm.y[s])    
+    JuMP.@constraint(dm.model,
+        [r=1:length(stages)],
+        dm.open[stages[r]] >=
+        sum(dm.y[s] for s in stages[r]) - length(stages[r]) + 1)
+end
+
+function ridetransitconstraint(dm::DesignModel)
+    routeoptions = Dict{Tuple{Int, Int}, Any}()
+    for u=1:dm.np.nstations, v=TransitNetworks.dests(dm.np, u)
+        # TODO: Replace with stages that are ranked higher than 
+        # no-transit option in preference ranking
+        routeoptions[u,v] = 
+            unique(TransitNetworks.expandedstages(dm.np,u,v)[2])
+    end 
+    JuMP.@constraint(dm.model,
+        [u=1:dm.np.nstations, 
+         v=TransitNetworks.dests(dm.np,u),
+         r=1:length(routeoptions[u,v])],
+        dm.ride[u,v] >= dm.open[routeoptions[u,v][r]])
+    JuMP.@constraint(dm.model,
+        [u=1:dm.np.nstations,
+         v=TransitNetworks.dests(dm.np,u)],
+        dm.ride[u,v] <= 
+        sum(dm.open[routeoptions[u,v][r]] 
+            for r in 1:length(routeoptions[u,v])))
+end 
+
+function budgetconstraint(dm::DesignModel)
+    costs = [sum([TransitNetworks.haversinedistance(dm.np, 
+                    seg[i], seg[i+1]) for i in 1:length(seg)-1])
+                for seg in dm.np.segments]
+    JuMP.@constraint(dm.model,
+        dot(costs, dm.y) <= dm.np.budget)
+end 
+
+function segmentmatchingconstraint(dm::DesignModel)
+    for s1 in 1:dm.np.nsegments, s2 in (s1+1):dm.np.nsegments
+        if (sort(dm.np.segments[s1]) == sort(dm.np.segments[s2]))  
+            JuMP.@constraint(dm.model, dm.y[s1] == dm.y[s2])
+        end
+    end 
+end 
diff --git a/src/model/designmodel.jl b/src/model/designmodel.jl
new file mode 100644
index 0000000..7634e07
--- /dev/null
+++ b/src/model/designmodel.jl
@@ -0,0 +1,42 @@
+mutable struct DesignModel
+    np::TransitNetworkProblem
+    model::JuMP.Model
+    
+    auxinfo::Dict{Symbol,Any}
+    results::Dict{Symbol,Any}
+
+    y::Vector{JuMP.Variable}
+    ride::JuMP.JuMPDict{JuMP.Variable}
+    open::JuMP.JuMPArray{JuMP.Variable,1,Tuple{Array{Array{Int64,1},1}}}
+
+    function DesignModel(
+        np::TransitNetworkProblem,
+        model::JuMP.Model
+    )
+        new(np, model, Dict(), Dict())
+    end
+end
+
+function DesignModel(
+        np::TransitNetworkProblem,
+        solver = Gurobi.GurobiSolver(OutputFlag=0),
+        warmstart = zeros(0,0)
+    )
+    dm = DesignModel(np, JuMP.Model(solver=solver))
+
+    # Variables
+    dm.y    = designvariable(dm)
+    dm.ride = ridershipvariable(dm)
+    dm.open = openroutevariable(dm)
+
+    # Constraints
+    routeopenconstraint(dm)
+    ridetransitconstraint(dm)
+    budgetconstraint(dm)
+    segmentmatchingconstraint(dm)
+
+    # Objective
+    ridershipobjective(dm)
+
+    dm
+end 
diff --git a/src/model/objective.jl b/src/model/objective.jl
new file mode 100644
index 0000000..ff80728
--- /dev/null
+++ b/src/model/objective.jl
@@ -0,0 +1,3 @@
+function ridershipobjective(dm::DesignModel)
+    JuMP.@objective(dm.model, Max, sum(dm.ride))
+end 
diff --git a/src/model/variables.jl b/src/model/variables.jl
new file mode 100644
index 0000000..4fb5758
--- /dev/null
+++ b/src/model/variables.jl
@@ -0,0 +1,20 @@
+function designvariable(dm::DesignModel)
+    JuMP.@variable(dm.model,
+        y[1:dm.np.nsegments], Bin)
+    y
+end
+
+function ridershipvariable(dm::DesignModel)
+    JuMP.@variable(dm.model, 
+        ride[u=1:dm.np.nstations, v=TransitNetworks.dests(dm.np,u)], Bin)
+    ride 
+end
+
+function openroutevariable(dm::DesignModel)
+    validod = vec([isassigned(dm.np.stage, u,v) 
+        for u in 1:dm.np.nstations, v in 1:dm.np.nstations])
+    stages = unique(vcat(dm.np.stage[validod]...))
+    JuMP.@variable(dm.model,
+        open[stages], Bin)
+    open
+end