sz = 50;
d = Array[Abs[#1 - #2] &, {sz, sz}];
vd = N@Rest@Flatten@d;
fs = Array[f, {sz, sz}];
vfs = Rest@Flatten@fs;
uv = ConstantArray[1, sz];
on = Normalize[RandomReal[1, sz], Total];
tw = Normalize[RandomReal[1, sz], Total];
eqs = {fs[[1, 1]], -on - fs.uv, -tw - uv.fs} /. Solve[-1. == fs /. List -> Plus, f[1, 1]][[1]] // Flatten;
{b, m} = CoefficientArrays[eqs, vfs];
AbsoluteTiming[
 out = LinearProgramming[vd, m, b, Method -> "InteriorPoint"];
 vd.out
]