What you are trying to do is very well possible. The problem you are encounting is, that your fitting algorithm crashes due to a singular matrix inversion. You can resolve that problem in a couple of ways. The easiest is to limit the amount of iteration to find the fitting curve. So this script:
a=1
b=1
FIT_MAXITER = 1
f(x) = a*x + b
ti = sprintf("%.2fx+%.2f", a, b)
fit f(x) 'data' using 1:2 via a, b
plot [0:3] f(x) t ti, "data" w l
should do exactly what you are aiming for.
Note that the singular matrix inversion problem should not arise when your data is noisy or your setup function does not have the exact structure as your data. For example this
f(x) = a*x**2 + b
function should work just fine without limiting the number of iterations.
Further ways to control the fitting process are described in the gnuplot documentation (gnuplot.pdf or help set fit).