From first principles, I know that the best reduced-order mathematical model that describes my complicated physical system is
$$\dot{x}(t) = \alpha f(u,t) + \beta u(t)$$
where $\left(\alpha,\beta\right)\in \mathbb{R}$ are known constants.
I have collected a bunch of step-response data from simulating my complicated physical model, by varying the magnitude of my step input $u(t)$. Note that there is no noise in the input or output data. From all the step-response plots, I deduce that my system behaves like a classic first-order system, i.e. the debiased response looks something like
$$x(t) - x(0) = x_\mathrm{final}(1 - e^{-t/\tau})\, u(t)$$
For this reason, I suspect that my (reduced order) model is an LTI system wherein $f$ from the model hypothesis is at best a weak function of $t$ and furthermore, has only a negligible coupling to $u(t)$ (maybe I am wrong?)
a) What is the best procedure for identifying the $f$ component in the RHS of the model? I am confused because this equation form does not (at-present) look like a standard state-space model. Is it expected to reformulate this into a standard LTI transfer function or state-space equation?
b) If the above option is not feasible, what would be the best system-ID procedure that can identify any other suitably reformulated model that can capture the model dynamics "sufficiently" enough, i.e. predict $(x(t) \text{ given } (\alpha,\beta, u(t))$?
If continuous system ID is not possible, I am happy to accept a discrete-time model. Please advise me of the procedure as I am a novice in system-ID and Dr. Ljung's textbook seems daunting to parse.
One can test whether a system might be nonlinear by applying a sinusoidal input signal, $u(t) = A_u\,\sin(\omega\,t + \phi_u)$. If the system is LTI then the output should be a transient plus also a sinusoidal signal, but with a different amplitude and phase $y(t) = A_y\,\sin(\omega\,t + \phi_y)$. If the system is stable, then this transient should die out over time. This seems to be the case according to your step responses. For LTI systems scaling this input by a constant should also result into the same scaling of the output, $u'(t) = \alpha\,u(t)\,\to\,y'(t) = \alpha\,y(t)$. For nonlinear systems these input/output relations often do not hold.
If the system seems to be LTI then it is possible to fit a transfer function onto this amplitude and phase information. Namely this information can be represented as a frequency response function.
Without further knowledge of the system it will be hard to try and fit a nonlinear model. However since you expect it to be or close to LTI, then fitting such system should, for a certain range of inputs, capture the system dynamics sufficiently.