Efficient formulation for binary integer linear programming
Integer linear programming

Let me suggest another way of formulating this with ILP that might be worth trying.  Define combination to mean a list of all of the balls in a single box.  For instance, the combination might be (7,15) meaning that the box contains ball 7 and ball 15.  Of course, we can enumerate all legal values for the combination, i.e., for the contents of a single box.  There will be at most $1+N_S+N_B+N_S(N_S-1)/2+N_B(N_B-1)/2+N_B N_S(N_S-1)/2$ different combinations (fewer in practice due to the constraints on the difference of weights and the total weight of a box).

Now introduce a binary variable $x_{ij}$ that is 1 if the $i$th box contains combination $j$.  Here $j$ is an index that ranges over all possible legal choices for the combination, i.e., for the contents of a single box.  (Don't include any illegal combinations.)

We get some linear inequalities from this:


For each box, we must select one combination for it to contain: $\sum_j x_{ij} \le 1$.
Each ball must be used exactly once: $\sum_i \sum_{j \in B_k} x_{ij}=1$, where $B_k$ is the set of legal combinations that include ball $k$.


Now you can test feasibility, i.e., test whether it is possible to pack all of the balls into some given number of boxes.

For example, suppose we have 5 big balls and 10 little balls and 10 boxes.  Then there will be at most about 300 legal combinations (probably fewer), so there will be about 3000 binary variables.  I don't know whether PuLP will be effective, but you could give it a try.

Exact cover

Even better yet, one can note that your problem is an instance of the exact cover problem.  The universe is the set of balls; each combination is a subset of the universe; and we want to find a way to cover the universe by a disjoint union of the combinations.

There are dedicated algorithms for the exact cover problem that might be faster than expressing this as an integer linear program and applying an ILP solver.