To impose the distance restriction, use a sparse index set of $(i,j)$ pairs rather than the full Cartesian product $I \times J$.
Also, you might consider omitting constraint $(2)$, which will naturally be satisfied unless the penalty for unmet demand is too small to encourage opening any facilities, and constraint $(5)$, which is logically implied by $(3)$ ...
Following is a possible way of its implementation:
First, you define a function that takes an OD matrix, solves the GUROBI model, and returns the optimal locations.
// This part will prepare the gurobi model and change the parameters related to OD matrix
// Solve the model and return locations
You can add a constraint that says the number of vulnerable people assigned to a facility is at least a specified fraction of the total number of vulnerable people (where the fraction is set to 1 if you want to ensure all vulnerable people are assigned, or something less than 1 indicating your tolerance for leaving vulnerable people unassigned). If you use a ...