diff --git a/src/LabProcesses.jl b/src/LabProcesses.jl
index 31f15627660295d7d4d4b54fb7008565033ea8bb..b64e23d7139d020b56bfdebaf79e3d92bfcdd355 100644
--- a/src/LabProcesses.jl
+++ b/src/LabProcesses.jl
@@ -27,6 +27,9 @@ sampletime(p::AbstractProcess)  = error("Function not implemented for $(typeof(p
 control(p::AbstractProcess, u)  = error("Function not implemented for $(typeof(p))")
 measure(p::AbstractProcess)     = error("Function not implemented for $(typeof(p))")
 
+initialize(p::AbstractProcess)  = error("Function not implemented for $(typeof(p))")
+finalize(p::AbstractProcess)    = error("Function not implemented for $(typeof(p))")
+
 
 # Interface implementation Ball And Beam ====================================================
 
@@ -59,4 +62,9 @@ const conversion = 65535/20
 io2num(x) = x/conversion -10            # Converts from io to float
 num2io(x) = round(Int32,x*100 + 2050)   # Converts from regular number to io
 
+initialize(p::BallAndBeam)  = ccall((:comedi_start, comedipath),Int32,(Int32,), 0)
+finalize(p::BallAndBeam)    = ccall((:comedi_stop, "../c/comedi_bridge.so"),Int32,(Int32,), 0)
+initialize(p::BallAndBeamSimulator)  = nothing
+finalize(p::BallAndBeamSimulator)    = nothing
+
 end # module
diff --git a/src/controllers.jl b/src/controllers.jl
new file mode 100644
index 0000000000000000000000000000000000000000..ec62de6359bac5279120bf6fbe800e24b69d0449
--- /dev/null
+++ b/src/controllers.jl
@@ -0,0 +1,42 @@
+"""
+	y,u,r = run_control(process, sysFB[, sysFF]; duration = 10, reference(t) = sign(sin(2π*t)))
+Perform control experiemnt where the feedback and feedforward controllers are given by
+`sysFB` and `sysFF`, both of type `StateSpace`. See [`fbdesign`](@ref) [`ffdesign`](@ref).
+
+`reference` is a reference generating function that accepts a scalar `t` (time in seconds) and outputs a scalar `r`, default is `reference(t) = sign(sin(2π*t))`.
+
+The outputs `y,u,r` are the beam angle, control signal and reference respectively.
+
+"""
+function run_control(P::AbstractProcess,sysFB, sysFF=nothing; duration = 10, reference = t->sign(sin(2π*t)))
+	h 		= sampletime(P)
+	y       = zeros(0:h:duration)
+	u       = zeros(0:h:duration)
+	r       = zeros(0:h:duration)
+
+	sysFB   = minreal(sysFB)
+	stateFB = init_sysfilter(sysFB)
+	if sysFF != nothing
+		sysFF   = minreal(sysFF)
+		stateFF = init_sysfilter(sysFF)
+	end
+
+	function control(i)
+		rf = sysFF == nothing ? r[i] : sysfilter!(stateFF, sysFF, r[i])
+		e  = rf-y[i]
+		ui = sysfilter!(stateFB, sysFB, e)
+	end
+
+	initialize(P)
+	for (i,t) = enumerate(0:h:duration)
+		@periodically h begin
+			y[i]       = measure(P)
+			r[i]       = reference(t)
+			u[i]       = control(i) # y,r must be updated before u
+			control(P, u[i])
+		end
+	end
+	control(P, 0)
+	finalize(P)
+	y,u,r
+end