Consider a Hammerstein nonlinear dynamical system of the form
$\mathbf{\dot{x}} = \mathbf{Ax} + \mathbf{Bu}$,
where the non-linearity is in the control term $\mathbf{u}$, and has a piecewise conditional form, viz.,
$\mathbf{u} = \left\{\begin{array}{cc} \mathbf{x}-\mathbf{\delta}, & x_i \ge \delta_i, \\ 0, & |x_i| < \delta_i \\ \mathbf{x} + \mathbf{\delta}, & x_i \le -\delta_i.\end{array}\right.$
and where the subscript $i$ represent's the $i^{\rm{th}}$ component of the vector (for some fixed $i$).
Such a system represents a symmetric dead-zone nonlinearity and appears often in engineering applications. Now let $\delta_i$ be a random variable with some distribution (on non-negative support), and consequently $\mathbf{x}(t)$ is a process parameterized by a random variable: $\mathbf{x}(t;\delta)$.
My question is this: what is the best way to model the uncertainty in this dynamical system? The random variable appears in the conditional term of the piecewise definition; therefore any treatment must consider the effect of uncertainty not just in the moment created by the $\mathbf{x}-\delta$ term, but also in the location of the breakpoint.
In the past, I have used non-intrusive polynomial chaos with some success to model this effect. However, I am curious if there is another possibly more clever and accurate technique.