Your problem is you have a multiple input multiple output system which couples the elements of your sweep, thus we indeed expect the edge effects of your zero-phase digital filter to work their way inwards. You can consider that at each sweep you get a vector of inputs, rather than a timeseries of data, with the transfer function of your system being: $$\mathbf{y}=\mathbf{x}-\left(\frac{1}{1-z^{-1}}\mathbf{I}\right)\mathbf{H}\left(z^{-1}\mathbf{I}\right)\mathbf{y}=\left(\mathbf{I}+\left(\frac{z^{-1}}{1-z^{-1}}\mathbf{I}\right)\mathbf{H}\right)^{-1}\mathbf{x}$$ Where lower case boldface are vectors, and upper case boldface are matrices ($\mathbf{I}$ being the identify matrix). That matrix $\mathbf{H}$ (your FIR filter) is coupling everything, which is why it's expected that the edge effects will work their way inwards. To assess stability, you'd have to solve: $$\mathbf{det}\left(\mathbf{I}+\left(\frac{z^{-1}}{1-z^{-1}}\mathbf{I}\right)\mathbf{H}\right)=0$$ This is much more complex than things need to be.
As a solution, redesign your control loop as per my figure below, which will guarantee stability and prevent the edge effects creeping inwards, with the transfer function given by: $$\mathbf{e}=\mathbf{H}\left(L[z]z^{-1}\mathbf{I}\right)(\mathbf{e}+(\mathbf{x}-\mathbf{e}+\mathbf{d}))=\mathbf{H}\left(L[z]z^{-1}\mathbf{I}\right)(\mathbf{x}+\mathbf{d})$$ $$\mathbf{y}=\mathbf{x}-\mathbf{e}+\mathbf{d}=\left(\mathbf{I}-\mathbf{H}\left(L[z]z^{-1}\mathbf{I}\right)\right)(\mathbf{x}+\mathbf{d})$$ Where $L$ is the transfer function of a low pass filter to provide some averaging of uncorrected sweeps (probably a moving average filter will do, but I leave that to you), and I have included a disturbance, $\mathbf{d}$, to be more rigorous because I presume your DAC and ADC with have some inaccuracy that is uncorrelated with the signal. This system clearly has no poles, except possibly in the low pass filter, $L$, so design a stable low pass filter, and you're guaranteed a stable system. Hopefully that final transfer function is quite intuitive in that, ignoring the disturbance, $\mathbf{d}$, it takes $\mathbf{x}$ and subtracts the component you don't want to give you $\mathbf{y}$, which is exactly what you want it to do. Note that the role of the low pass filter is just to give some averaging between sweeps to attenuate their noise, which I guess you are aware of that you need insofar that one sweep of data will not provide an accurate correction, which is why I presume you used that integrator. On the first iteration, I would initialise the low pass filter so that its output immediately is equal to the first sweep signal, like what your integrator did.
