Grey-box system ID for a set of DAEs

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 grey-box system-id approach that shall work well to arbitrary excitations in all future simulations. (NOT merely a curve-fitting procedure to match a specific excitation input-output dataset)