# Procedure to identify characteristic properties of unknown functions in a DAE model

I have a system of 1st order odes given by

$$\dot{x_1}(t) = \alpha_1 f_1(x_1,t) + \beta_1 u(t) \\ \dot{x_2}(t) = \alpha_2 f_2(x_2,t) + \beta_2 u(t)$$

They are constrained by an algebraic equation

$$x_1(t) + x_2(t) = k$$

where $\left( \alpha_1,\alpha_2, \beta_1,\beta_2 , k \right) \in \mathbb{R}$ are known constants (i.e. parameters). $f_1(t)$ and $f_2(t)$ are both unknown.

Starting from a rich set of input-output noise-free data available from simulating a complex proxy system, what would be the best procedure to identify (even a subset of repeatable/characteristic properties) the unknown possibly time-varying functions $f_1(x_1,t)$ and $f_2(x_2,t)?$ I am almost certain that $f_1(x_1,t)$ and $f_2(x_2,t)$ are both linear.

I am looking for a suitable approach that shall work well to arbitrary excitations in all future simulations. (NOT merely a curve-fitting procedure to match a specific input-output dataset)