I think you are almost in the right way. What you need to do is not just replace the variance on the diagonal of cov matrix, but use multidimensional GARCH model to predict the whole covariance matrix. Such as BEKK model.
Portfolio optimization has nothing to do with historical data or forecasted data. If you have a return vector and a covariance matrix, the optimization can be done. So if you can predict the covariance matrix and return vector by multidimensional GARCH model, everything is done.
I recommend a paper by Woodridge. Perhaps it is not so close but it may be helpful. At least it gives some basic idea. 
Hope Help~