The fundamental problem is that you have an expression raised to a negative power, and if that expression can ever be 0, that gives you an infinite contribution. Your expression includes (a-x) so if your a parameter can ever exactly equal any of the x values, you generate the infinity. If (a-x) is negative (so if a is less than any x value) then you get a complex result, since you would then have negative to a fraction.
The reason the problem showed up on the initial call, is that your initial guess of [1 1] had a value that exactly matches the upper bound of 1 for the first parameter. The internal processing of bounds pushed the guess to an interior point, and the interior point happened to be 0.55 ... which happens to be exactly one of the x values, so you got 0^negative giving you the infinity.
In the below, I set the lower bound for the first parameter to be the maximum x value; that could potentially still give you 0^ but since it avoids touching the bound itself it moves the guess to a valid range.
fy = @(a,b,x) 0.001*((a-x)./(a.*(1-x))).^(-b.*a./(1-a));
x = [0,0.1, 0.2, 0.3, 0.4, 0.46, 0.55, 0.6];
y = [0.001, 0.00111499, 0.0011926, 0.0013699, 0.00161633, 0.00192075, 0.00274991, 0.00357156];
fun = @(c) y - fy(c(1),c(2),x)
residue = @(c) norm(y - fun(c));
lb = [0 0];
ub = [1 90];
lb(1) = max(x);
guess = [1;1];
B = fmincon(residue, guess, [],[],[],[], lb, ub);
fprintf('PhiMax = %.15f\nVisco = %.15f\n', B)
xv = linspace(min(x), max(x));
figure(1)
plot(x, y, 'pg')
hold on
plot(xv, fy(B(1),B(2),xv), '-r')
hold off
grid
