There does not seem to be a standard way to deal with missing data in the context of the exponential smoothing family of models. In particular, the R implementation called ets in the forecast package seems to just take the longest subsequence without missing data, and the book "Forecasting with Exponential Smoothing" by Hyndman et al. does not seem to talk about missing data at all.
I would like to do a little more, if my users explicitly ask me to (and if the missing data don't occur too close together or in too many periods that are exactly one season apart). In particular, what I have in mind is the following. During simulation, whenever I would encounter a missing value $y_t$, I would substitute the current point forecast $\tilde y_t$ for $y_t$, so that $\varepsilon_t = 0$. This would, for example, make the data point not be considered for the parameter optimization process.
Once I have a reasonable fit for the parameters, I can estimate the standard deviation of the errors (assumed to be normal with mean $0$) and verify that using values for $\epsilon_t$ generated from that distribution do not decrease the likelihood by a big factor. I would use such values for forecasting (using simulation) as well.
Are there known pitfalls with this method?