Let's get some potential obstacles out of the way first. One theme that's central to all renormalization narratives is to relate the bare theory to physical observables, which is often manifested as the study of the mathematical relations between bare and renormalized couplings. This is what was studied from the beginning of QFT renormalization program, and this central point hasn't been changed ever since, in particular renormalization group perspective doesn't change it. What was troublesome at the time was the insistence on sending cutoff to infinity, which is more than often impossible(e.g. Landau pole), but today people no longer insist on that. Instead, we keep the weaker requirement that cutoff may stay but the low energy physics must be insensitive to the cutoff, for example we can construct our QFT such that S-matrix elements depend on $O(\Lambda^{-1})$, so we might as well just take $\Lambda\to\infty$, but beware this by no means suggest the cutoff can be consistently removed from all parts of our QFT(again, e.g. Landau pole). When the cutoffs can genuinely be removed, we say we have a completion, this is the case for Yang-Mils theory, for example.

When we say "Wilsonian renormalization group (philosophy)", especially when used in the sense of "We only understand what QFT is after Wilson"(somewhat an exaggeration in my opinion), we are typically referring to the realization that renormalization is a sort of coarse graining process (and indeed technically very close to what OP described in the 2nd last paragraph), that a completion of QFT corresponds to the existence of fixed points(modulo limit cycles that are usually not of interest), that nonrenormalizable operators naturally die out at low energy, and that renormalization is intimately connected to critical phenomena etc. These equip us with not only much better intuitions but also very powerful new tools for calculations.

The Wilson(-Polchinski) integrating-momentum-shell approach you referred to in your main post is just one of the many ways to do coarse graining, it's not absolutely necessary to use this approach to talk about Wilsonian philosophy, but it is a convenient one.

All said above involve changing the cutoffs, and are more or less attacking the same problem, with different levels of sophistication: **find a family of theories/Lagrangians labeled by cutoffs, such that they give precisely the same prescribed physics at some energy scale lower than the cutoff. **Then by changing the cutoff, your coupling constants, as functions of cutoff, will form a trajectory. If you can make your cutoff arbitrarily large without any inconsistency, congratulations, you have a continuum completion; if not, it's not necessarily a catastrophe, maybe you are dealing with a situation with a physical cutoff, such a solid state lattice system, or that at certain scale a qualitatively different unknown theory has to step in.

As for Callan-Symanzik equations, there are several Callan-Symanzik-like equations, and I suspect all of them have been referred to as "Callan-Symanzik equation" in literatures.

For example, one may ask, as discussed at length above, if we keep the renormalized couplings(at a fixed energy scale) fixed, how would bare couplings change as a function of cutoff. One may also ask if we keep the bare couplings constant, how would renormalized couplings(at a fixed energy scale) change as a function of cutoff. Clearly these two are mathematically equivalent, they both yield Callan-Symanzik-like equations, and both are in accordance with the philosophy described in the lengthy discussion above.

The last sort keeps bare couplings fixed, cutoff fixed, but ask how the renormalized couplings change as the reference energy scale changes, this is what you are looking for, but unfortunately I myself am not sure if this is equivalent to the former two Callan-Symanzik-like equations, all I know is that the beta functions of all these three equations have the same power series expansion coefficients for the two leading terms, so at least there should be some weak equivalence in the practical sense. I'll leave it to who knows better.

-------------------------------------------------------------------------------------------------------------

Update:

1. I should add that in the famous review by Wilson and Kogut on RG and epsilon expansion, they call the curve that is produced by taking $\Lambda$ to infinity "canonical curve", while they call the curve produced by changing the coarse-graining/measurementl scale ($\Lambda_N$ or $\mu$ in OP's noatation) "RG trajectory", so they made it quite explicit that the two are in principle two different things.

2. I now understand exactly why perturbatively the canonical curve and RG trajectory have the same form of $\beta$ function up to the first two orders, it's very lucidly explained in the brilliantly written book "Quarks, gluons and lattices" by Michael Creutz:

For a theory with only a dimensionless coupling, any dimensionless observable H must have the functional dependence on cutoff scale $\Lambda$, measurement scale $\mu$ and the bare coupling $g_0(\Lambda)$ in the particular form of $H(\mu,\Lambda, g_0)= H(\frac{\mu}{\Lambda},g_0)$, by simple dimensional analysis. This immediately gives

$$\Lambda \frac{\partial}{\partial\Lambda} H=-\mu\frac{\partial}{\partial\mu}H$$

The renormalized coupling is always defined as the value of some observable at certain scales, so it also satisfies

$$\Lambda \frac{\partial}{\partial\Lambda} g_R=-\mu\frac{\partial}{\partial\mu} g_R$$

On the other hand, as we push the cuttoff to infinity, the $g_R$ should gradually become a constant with respect to $\Lambda$. (Of course, if keeping $g_R$ fixed at a certain scale is one your renormalization conditions, then this is automatically true for all values of $\Lambda$.) This means

$$0=\Lambda \frac{d}{d\Lambda} g_R=\Lambda \frac{\partial}{\partial\Lambda} g_R+\frac{\partial g_R}{\partial g_0} \Lambda \frac{d}{d\Lambda} g_0.$$

Put together them we have

$$\beta^R(g_R)\equiv \mu\frac{\partial}{\partial\mu} g_R=\frac{\partial g_R}{\partial g_0} \Lambda \frac{d}{d\Lambda} g_0 \equiv\frac{\partial g_R}{\partial g_0} \beta(g_0).$$

On top of all these, to deserve the name "renormalized **coupling**", $g_R$ should be defined such that it equals $g_0$ at leading order, so $g_R= g_0+bg_0^3+O(g_0^5)$.

At this stage one can apply the standard calculation to show why the first two coeffiencients in the expansion $\beta^R(g_R)$ (in $g_R$) are the same as those of $\beta(g_0)$ (in $g_0$).