You are correct that Cauchy-Kovalevskaya does not apply directly to this problem, but there are other theorems that give sufficient conditions, provided that you make certain basic regularity assumptions.
For example, in the smooth involutive case, i.e., when the $Z_i$ (as well as the $\lambda_i$) are also sufficiently differentiable and satisfy the involutivity condition that $[Z_i,Z_j] = c^k_{ij}\,Z_k$ (summation convention assumed) for some smooth functions $c^k_{ij}=-c^k_{ji}$, then an obvious necessary condition on the $\lambda_i$ for local solvability of the system $Z_if = \lambda_i$ is
$$
Z_i(\lambda_j)-Z_j(\lambda_i) = c^k_{ij}\lambda_k\,\qquad \forall i,j
\quad\,1\le i,j\le r,\tag1
$$
since both sides would equal $[Z_i,Z_j]f$ for any (local) solution $f$. It's a consequence of the Pfaff-Darboux Theorem (appropriately interpreted) that, if condition (1) is satisfied, then, every $p\in U$ has an open neighborhood $V\subset U$ on which there exists an $f$ satisfying the given system $Z_if = \lambda_i$. However, because the system is linear inhomogeneous, there is a simpler proof of this result that only uses the Frobenius Theorem. I give that proof in a remark at the end of this answer.
This can be extended to more general situations: Suppose that the distribution $D$ spanned by the $Z_i$ is not involutive, but that we can choose new $Z_a$ for $r<a\le r'$ such that $Z_1,\ldots, Z_{r'}$ are linearly independent and give a basis for the sections of the distribution $D'$ of constant rank $r'$ spanned by the $Z_i$ and $[Z_i,Z_j]$ for $1\le i,j\le r$. Then writing $[Z_i,Z_j] = c^a_{ij} Z_a$ (where now, the index $a$ runs from $1$ to $r'$), we see that a necessary condition for solvability is that there exist functions $\lambda_a$ for $r<a\le r'$ (necessarily unique) that satisfy
$$
Z_i(\lambda_j)-Z_j(\lambda_i) = c^a_{ij}\lambda_a\,\qquad \forall i,j\qquad\,1\le i,j\le r,\tag2
$$
(where the implied summation on $a$ runs from $a=1$ to $a=r'$)
and that we would have to have $Z_af = \lambda_a$ for $1\le a\le r'$.
This is now a bigger system (the 'prolongation' of the original system), and we can check whether $D'$ is involutive. If not, we expand (aka,
'prolong') the system again. If we can repeat this prolongation process until we reach a system $D''$ for some $r''>r$ that is involutive, then we can apply the first criterion to determine local solvability of the original system.
Along the way, we might run into situations where the new $\lambda_a$ one needs to find don't exist on some open set, in which case, there won't be a solution to the original system. More complicated things can happen if the 'prolonged' distributions don't have constant rank, but, for most practical purposes, the above process will give an effective test as to whether and where in $U$ (local) solutions $f$ exist.
Remark: I can sketch the proof of the basic existence theorem mentioned above, based on the Frobenius theorem.
First, suppose that $[Z_i,Z_j] = 0$ for $1\le i,j,\le r$. Then, by the Simultaneous Flow Box Theorem, for each $p\in U$ there exists an open $p$-neighborhood $V\subset U$ on which there exists a $p$-centered coordinate chart $z = (z^1,\ldots,z^n)$ such that $Z_i = \partial/\partial z^i$ on $V$ for $1\le i\le r$. On $V$, the given system becomes $\partial f/\partial z^i = \lambda_i$, and the condition (1) becomes $\partial\lambda_i/\partial z^j = \partial\lambda_j/\partial z^i$, and, by the usual integration formula, it is clear that, when this holds, there will be a smaller $p$-neighborhood $V'\subset V$, on which there exists a function $f$ satisfying $\partial f/\partial z^i = \lambda_i$ for $1\le i\le r$.
Now consider the more general case in which $[Z_i,Z_j] = c^k_{ij}\,Z_k$. By the Frobenius Theorem, every point $p\in U$ will have a $p$-neighborhood $V$ with a $p$-centered coordinate system $z = (z^1,\ldots,z^n)$ such that $Z_i = a_i^k\,\partial/\partial z^k$ on $V$ for some smooth functions $a_i^k$ on $V$. Moreover, the $r$-by-$r$ matrix $a = (a_i^k)$ is invertible because the $Z_i$ for $1\le i\le r$ are linearly indepdendent on $U$. Consequently, we can write $\lambda_i = a_i^k\mu_k$ for some functions $\mu_k$ ($1\le k\le r$) on $V$.
Now, we can expand $[Z_i,Z_j] = c^k_{ij}\,Z_k$ to get a formula for $c^k_{ij}$ in terms of the $a_i^k$ and their partials with respect to the $\partial/\partial z^j$. Using this formula and some index juggling, we see that
$$
Z_i(\lambda_j)-Z_j(\lambda_i) - c^k_{ij}\lambda_k
= a_i^ka_j^l\left(\frac{\partial \mu_k }{\partial z^l}-\frac{\partial \mu_l}{\partial z^k}\right).
$$
Since the given system $Z_if=\lambda_i$ is equivalent to $\partial f/\partial z^k = \mu_k$, we see that we are reduced to the case $c^k_{ij} = 0$, which has already been treated.
