Fitting an unknown cost function in python via polynomial models, best way for cross validation

I do fit an ODE to some data (coupled diffusion and adsorption). However, I do not know the physical relationship between the parameters.Unfortunately, one parameter changes during integration (makes it hard to use constraints). Therefore, I do use a bounded polynomial model in order to calculate permeability similar to:

def get_diff_coeff_13(p, conc_salt, ph, ratio):
    """Calculates the lumped permeability from general logistic model.
       Permeability is limited between a (lower) and k (upper)
    :param p: parameter to optimize over
    :param conc_salt:
    :param ph:
    :param ratio:
    :return:
    """
    x = p[0]*conc_salt + p[1]*ph + p[2]*ratio + p[3]*ph*conc_salt + p[4]*ratio*ph
    a = 0 # lower bound
    k = 1e-9 # upper bound
    c = 1
    b = 0.01 # determines the slope of the linear region
    m = 1
    n = 1
    diff = a + (k-a)/(c+np.exp(-b*(x-m)))**(1/n)
    return diff

As you might imagine from the name of the function. I do have quite a lot of combinations of polynomials. Currently my best idea is to use a brute force approach and fit all combinations of my polynomial within a loop and log residual sum of squares.

How would you approach such a problem ? Are there more elegant methods than using a brute force approach ?