I think my question can be posed in many PDE that admit stationary solutions or travelling waves. I will simply describe the setting that I am more familiar with.

Consider the viscous Burgers' equation in $\mathbb R$ with a forcing term: $$ v_t-v_{xx}=(v^2)_x +f_x,\qquad (x,t)\in \mathbb R\times [0,\infty) $$ $$ v(x,0)=v_0(x), \qquad x\in\mathbb R,$$ with kink-like boundary conditions at infinity: $v_0-\tanh(x)\in L^2(\mathbb R)$. The function $\tanh(x)$ is special because it is a stationary solution of the equation, together with all its translations. The only think to know is that $\tanh(x)\rightarrow \pm 1$ exponentially as $x\to\pm\infty$, and $\partial_x \tanh(x)=\operatorname{sech}^2(x)$. The forcing term $f$ lies in $L^2(\mathbb R\times (0,\infty)).$

The equation is globally well-posed in the affine space $L^2(\mathbb R)+\tanh$ by standard arguments.

Now, assume that the forcing term $f$ is small enough in $L^2(\mathbb R\times (0,\infty))$, and that initial datum is a small perturbation of some translation of $\tanh$, i.e., $|v_0|_{L^2}$ is small, where $$ |z|_{L^2}:=\min_{\gamma\in \mathbb R}\|z-\tanh(x-\gamma)\|_{L^2}. $$ Then, it is possible to prove that the quantity $|v(t)|_{L^2}$ is small for all $t>0$, and there exists a function $\alpha:[0,\infty)\mapsto \mathbb R$ such that $v=w+\tanh(x-\alpha(t))$ with the estimates $$ \|w\|_{L^\infty_tL^2_x}+\|w_x\|_{L^2_tL^2_x}\lesssim \|w_0\|_{L^2}+\|f\|_{L^2_tL^2_x}.\qquad (1) $$

Proof (sketch). To show that, the basic argument uses some modulational analysis: call $w(x,t):=v(x,t)-\tanh(x-\alpha(t))$, where we fix $\alpha(t)$ using the orthogonality condition $$ \int_{\mathbb R}\operatorname{sech}^2(x-\alpha(t))w(x,t)\,dx=0 \quad\forall t\geq 0.\qquad\qquad(2) $$ In this case, the function $w$ satisfies $$ w_t-w_{xx}-2(\tanh(x-\alpha)w)_x=(w^2)_x+f_x+\alpha'\operatorname{sech}^2(x-\alpha).\qquad (3) $$ Multiplying by $w$ and integrating in space (as for the usual energy estimates) we obtain $$ \frac 12\frac d{dt}\|w(t)\|_{L^2}+Q_\alpha(w)\leq\|f(t)\|_{L^2_x}\|w_x\|_{L^2}, $$ where the last term vanishes thanks to the orthogonality condition $(2)$, the term $(w^2)_x$ vanishes when tested against $w$, and $Q$ is the quadratic form $$ Q_\alpha(w):=\int_{\mathbb R}w_x^2-\operatorname{sech}^2(x-\alpha)w^2\,dx. $$ Moreover (there is some work to do to show this) one can prove that, if $w$ satisfies the orthogonality condition $(2)$, then $$ Q(w)\geq 0. $$ Thus, integrating in time, we get the nice estimate $$ \|w\|_{L^\infty_tL^2_x}+\|w_x\|_{L^2_tL^2_x}\lesssim \|w_0\|_{L^2}+\|f\|_{L^2_tL^2_x}. $$

Now, if you made it alive until here, thank you a lot! My problem now is:

Question: How to prove estimate $(1)$ for arbitrary, large data $v_0$ and $f$?

The main problem is that $\alpha(t)$ is not uniquely determined by the orthogonality condition $(2)$ if the solution is large in the $|\cdot|_{L^2}$ norm. This leads to the fact, for instance, that the function $\alpha$ a priori is not necessarily continuous. If I knew that $\alpha$ satisfying $(2)$ exists and is continuous, then I think I can prove the statement in a similar way by approximation, but I don't know how to prove that such a function $\alpha$ is continuous. In fact, if I pick $\alpha$ arbitrarily, it is perfectly possible to have something non-continuous. In the small data case, the continuity of $\alpha$ follows from equations $(2)$ and $(3)$:

Sketch. Differentiating $(2)$ in $t$ and substituting $w_t$ using $(3)$, one obtains: $$ |\alpha'(t)| \left |\int_{\mathbb R}\operatorname{sech}^2(x-\alpha)[\operatorname{sech}^2-w_x]dx\right|\leq\dots, $$ where the right hand side has good terms.

I feel like this is typical of modulational analysis for large data, but I have no idea how to treat this situation.