EDIT: Merged the two answers.
Here's the problem statement:
The input is $(V, \mathbf b, \ell)$, where $V=\{x_1,x_2,\ldots,x_n\}$ with each $x_i\in\{0,1\}^d$ (where $d$ is constant), vector $\textbf{b}=(b_1,b_2,\ldots,b_d)\in \mathbb{Z}_+^d$, and $\ell$ is an integer.
The problem is to decide whether $V$ can be partitioned into $\ell$ parts $V_1,V_2,\ldots,V_\ell$ such that $\forall~j,\sum_{x_i\in V_j} x_i \leq \textbf{b}$.
Lemma 1. The problem (with $d=O(1)$) reduces to a sparse language, so is not NP-hard unless P=NP.
Lemma 2. For $d=O(1)$ and $\ell=O(1)$, the problem can be solved in polynomial time via dynamic programming (EDIT: or an ILP of constant dimension).
Lemma 3. For $d=O(1)$ and $\max_j b_j=O(1)$, the problem can be solved in linear time (in the real RAM model) via an ILP of constant dimension.
This leaves open the complexity of the case when $\ell$ and $\max_j b_j$ are both allowed to grow (while $d=O(1)$). This case could be in P, or not in P but not NP-hard.
Proof of Lemma 1. One can faithfully re-encode any problem instance given in the form above into a representation that uses only $O(\log n)$ bits, as follows. For each of the $2^d$ possible vectors in $\{0,1\}^d$, write down the number of times that vector occurs in $V$. Write down $\ell$ and each entry of $\mathbf b$.
There are $2^d + d + 1$ numbers written, and each number is at most $n$, so the encoding size is $\Theta(2^d\log n)$ bits.) It follows that (for fixed $d$) the number of distinct instances of size $n$ is polynomial in $n$. This implies (details next) that the decision problem cannot be NP-hard unless P=NP.
In more detail, let $L$ denote the language for the problem encoded as first defined. Assume that $L$ is NP-hard. Define $L' = \{(f(V, \mathbf b, \ell), 1^{|V|}) : (V, \mathbf b, \ell) \in L\}$, where $f(V, \mathbf b, \ell)$ is the re-encoding of the instance as described above (into $O(\log |V|)$ bits), and $1^{|V|}$ is $|V|$ in unary.
$L$ reduces in poly-time to $L'$, so $L'$ is also NP-hard. But $L'$ is sparse (the number of elements of size $n$ is polynomial in $n$), so P=NP (by Mahaney's theorem). $~~~\Box$
Proof of Lemma 2.
The algorithm will solve the following more general problem.
The input is $(V, \mathbf B)$ where $V=\{x_1,\ldots,x_n\}$ (as before) and $\mathbf B = (B_1, B_2, \ldots, B_\ell)$ is a collection of $\ell$ vectors,
where each $B_j$ is in $\{0,1,\ldots,n\}^d$.
The problem is to decide whether there is a partition of $V$ into parts
$(V_1, V_2, \ldots, V_\ell)$ such that $\forall j,~\sum_{x_i\in V_j} x_i = B_j$.
This differs from the original problem in two ways. First, each part $V_j$ in the partition can have a different required sum $B_j$ (instead of $\mathbf b$). Second, equality must hold (the sum must equal $B_j$, rather being $\le \mathbf b$).
A given instance $(V, \mathbf b, \ell)$ of the original problem poly-time reduces to this one by considering all instances of this one of the form $(V, \mathbf B)$ where each $B_i$ is in $\{0,1,\ldots,n\}^d$ and $B_i \le b$. There are $O(n^{d\ell})$ such instances, which is polynomial since $d$ and $\ell$ are constant. The original instance will have a solution if and only if at least one of the latter instances does.
The algorithm solves the variant using the following dynamic program. Fix an instance $(V,\mathbf B)$ of the problem. Let $V=\{x_1,\ldots,x_n\}$ and $\mathbf B = (B_1,\ldots,B_\ell)$. So $x_i\in\{0,1\}^d$ and $B_j\in \{0,1,\ldots,n\}^d$ for each $i$ and $j$.
Define the signature of the instance $(V, \mathbf B)$ (following the proof of Lemma 1) to be obtained by replacing $V$ by the natural encoding of the function $f_V$ such that, for each $x\in\{0,1\}^d$, $f_V(x)$ is the number of times $x$ occurs in $V$. Note that any two instances with the same signature are equivalent, in that they have the same answer.
The signature is determined by $2^d + d\ell$ numbers, each in $\{0,1,\ldots,n\}$, so the number of possible signatures (among instances with the given parameters $d,\ell, n$) is $O(n^{2^d+d\ell})$, which (given that $d$ and $\ell$ are constant) is polynomial in $n$. Furthermore, the number of possible signatures of all "smaller" instances (meaning those where $n$ is smaller than the given $n$) is at most a factor of $n$ larger, so also polynomial in $n$.
Now we can proceed via dynamic programming. Fix the given instance $(V, \mathbf B)$. For the base case, if $V=\emptyset$ (that is, $n=0$), then the instance is feasible iff each $\mathbf B \equiv \mathbf 0$. Otherwise $x_n$ can go into one of the parts $V_j$ for $j\in\{1,2,\ldots,\ell\}$, such that removing $x_n$ from that part gives a solution to the problem obtained by subtracting $x_n$ from $\mathbf B_j$ (as long as $x_n\le \mathbf B_j$; if not $x_n$ cannot go in $V_j$).
For each of the (at most) $\ell$ possible ways of placing $x_n$, recursively determine whether the resulting subproblem has a solution. The original problem will have a solution if and only if one of these subproblems does.
Furthermore, memoize -- keep a cache of the answers to already-solved subproblems, indexed by signature. When encountering a given subproblem, check whether any subproblem with the same signature has already been answered, and if so use that answer instead of recursing. Because there are only polynomially many distinct signatures of smaller subproblems, the dynamic-programming algorithm will solve only polynomially many subproblems, so will run in time polynomial in $n$. $~~\Box$.
EDIT: As pointed out in the comments, an easier approach is (similarly to the proof of Lemma 3 below) to formulate an ILP of constant dimension, with a variable $y_{xj}$ for each $x\in\{0,1\}^d$ and $j\in\{1,\ldots,\ell\}$ representing the number of vectors in $V$ that are equal to $x$ and are assigned to part $V_j$.
The proof above raises the question of whether the variant defined in the proof is NP-hard when $\ell$ is not constant. (It might be, even if the original problem is not.)
Proof of Lemma 3.
The proof is an elaboration of @Gamow's answer in the comments.
Use the following equivalent problem formulation: Given a collection $V=(S_1, S_2, \ldots, S_n)$ of subsets of $\{1,2,\ldots, d\}$, where $d$ is constant, partition $V$ into a minimum number of parts, such that within each part each element $j$ occurs in more than $b_j$ subsets. Here is the algorithm.
Fix an input $V=(S_1, \ldots, S_n)$ of subsets of $\{1,2,\ldots, d\}$.
Let $\mathcal S_d$ denote the set of non-empty subsets of $\{1,2,\ldots,d\}$.
Let $\mathcal P_d$ denote the set of possible parts, that is, multi-subsets of $\mathcal S_d$ in which each element $j$ occurs at most $b_j$ times.
The goal is to partition the input $V$ into a minimum number of parts, each of which is in $\mathcal P_d$.
Note that $|\mathcal S_d| < 2^d$ and $|\mathcal P_d| \ll (bd)^{bd}$, so (as $b$ and $d$ are constant) there are constantly many possible subsets, and constantly many possible parts.
For each subset $S\in\mathcal S_d$ and possible part $p\in\mathcal P_d$, let $n(S, p)$ (either 0 or 1) denote the number of times that $S$ occurs in $p$.
For each possible subset $S\in\mathcal S_d$, count the number of times that $S$ occurs in the input $V$. Let $n(S, V)$ denote this number.
Finally, construct and solve the following integer linear program (ILP), with an integer-valued variable $x_p$ for each possible part $p\in\mathcal P_d$:
$$\text{minimize} \sum_{p\in\mathcal P_d} x_p \text{ subject to }$$
$$\begin{align}
(\forall S\in \mathcal S_d) && \sum_{p\in\mathcal P_d} n(S, p)\,x_p &{} = n(S, V) \\
(\forall p\in \mathcal P_d) && x_p &{} \in \{0,1,2,\ldots\} \\
\end{align}$$
For $\max_j b_j = O(1)$ and $d=O(1)$, the ILP has constant dimension (and a constant number of constraints) so can be solved in constant time in the real RAM model (by Megiddo's algorithm or subsequent improvements). The optimal value is the answer for the given instance $V$. $~~\Box$
The real RAM model does not account for bit complexity of arithmetic operations. But in each feasible solution to this ILP each coordinate has $O(\log n)$ bits, so presumably the ILP can be solved in polylog$(n)$ time in the standard RAM model,
so this approach should also yield a linear-time algorithm in the standard RAM model.
