  We analyse the far-infrared properties of $\sim$ 5,000 star-forming galaxies
at $z<4.5$, drawn from the deepest, super-deblended catalogues in the GOODS-N
and COSMOS fields. We develop a novel panchromatic SED fitting algorithm,
$\texttt{Stardust}$, that models the emission from stars, AGN, and infrared
dust emission, without relying on energy balance assumptions. Our code provides
robust estimates of the UV-optical and FIR physical parameters, such as the
stellar mass ($M_*$), dust mass ($M_{\rm dust}$), infrared luminosities
($L_{\rm IR}$) arising from AGN and star formation activity, and the average
intensity of the interstellar radiation field ($\langle U \rangle$). Through a
set of simulations we quantify the completeness of our data in terms of $M_{\rm
dust}$, $L_{\rm IR}$ and $\langle U \rangle$, and subsequently characterise the
distribution and evolution of these parameters with redshift. We focus on the
dust-to-stellar mass ratio ($f_{\rm dust}$), which we parametrise as a function
of cosmic age, stellar mass, and specific star formation rate. The $f_{\rm
dust}$ is found to increase by a factor of 10 from $z=0$ to $z=2$ and appears
to remain flat at higher$-z$, mirroring the evolution of the gas fraction. We
also find a growing fraction of warm to cold dust with increasing distance from
the main sequence, indicative of more intense interstellar radiation fields,
higher star formation efficiencies and more compact star forming regions for
starburst galaxies. Finally, we construct the dust mass functions (DMF) of
star-forming galaxies up to $z=1$ by transforming the stellar mass function to
DMF through the scaling relations derived here. The evolution of $f_{\rm dust}$
and the recovered DMFs are in good agreement with the theoretical predictions
of the Horizon-AGN and IllustrisTNG simulations.
