Support for fitting models with units to data with units

The astropy.modeling package currently contains a number of pre-defined models as well as ways for users to define their own models, and also contains a number of different fitting algorithms. However, neither the models nor the fitters currently work with values and arrays that have units attached to them (using astropy.units). It would therefore be nice to add support for units/quantities in astropy.modeling.

An example use case would be that users might want for example to fit a blackbody function (such as the one provided in astropy.analytic_functions, which is unit-aware) to fluxes in e.g. erg/cm^2/s/Hz or in mJy, and have the unit conversions be dealt with automatically.

There is already a prototype/work in progress implementation to add support for units to models, which can be found in astropy/astropy#3852, but this is not complete and does not include support for the fitting part, only the model evaluation.