r is a scalar?
How about something like:
model = direct_model(HiGHS.Optimizer())
set_silent(model)
@variable(model, x[i = 1:N] >= lb[i])
@variable(model, y[i = 1:N])
@objective(model, Min, l' * x)
@constraint(model, con[i in 1:N], 1 * x[i] - y[i] == 0)
@constraint(model, B * x - A * y .>= d)
for r in range(0, 1, length = n)
set_normalized_coefficient(con, x, fill(1 + r, N))
optimize!(model)
end
Now you only need to modify N coefficients in each step.