This is a really great question. I don’t know the answer yet. I struggle with this problem
every day when I design things in OpenSCAD, because OpenSCAD doesn't give your parametric model access to the results of the modeling operations in any way (!!)
I think that Antimony kind of suffers from the same problem, but because it's using distance fields rather than plain CSG, it does have a “blend” operation (https://github.com/mkeeter/antimony/blob/develop/py/fab/shap...) that attempts to fillet the union of two shapes. Now, Antimony has inversion (complement, more or less), and set A \ B is just A ∩ B', which by De Morgan's Laws is equivalent to (A' ∪ B)'. So what happens if you attempt to approximate a filleted A \ B with invert(blend(invert(A), B))? I don’t know. But damn, I’m sure tempted to try it!
Of course, even if that works, that won’t fillet a specific edge resulting from the subtraction. It will fillet all of them. Invariably, that is what I want.