In general, there cannot be a closed-form solution of a random coefficients VG model. The reason is the drift-restriction that needs to be imposed to ensure that the discounted price process is a martingale under the risk-neutral measure. Using the bank account as numeraire, the restriction is
$$ \frac{1}{\beta} > \theta + \frac{\sigma^2}{2} $$
where $\beta$ is the variance rate of the gamma subordinator, and $\theta$ and $\sigma$ are the drift and diffusion coefficient of the driving Brownian motion.
What can be done is to make the volatility of the driving process stochastic using a second (independent) process. So you could have a square-root process driving the volatility while the main process is subordinated to a gamma time. The book by Barndorff-Nielsen and Shiryaev, Chapter 12, contains some hints how to do that. The procedure would be analogous to the derivation of the regular VG characteristic function, but you would get the conditional characteristic function of the stochastic volatility model chosen instead of the normal distribution characteristic function from a regular Brownian motion. If you carefully select the volatility model, you can take the expectation with respect to the gamme time and get a closed-form solution of the characteristic function.