In this paper, a model is proposed for the optimal operation of multi-energy microgrids (MEMGs) in the presence of solar photovoltaics (PV), heterogeneous energy storage (HES) and integrated demand response (IDR), considering technical and economic ties among the resources. Uncertainty of solar power as well as the flexibility of electrical, cooling and heat load demand are taken into account. A p-efficient point method is applied to compute PV power at different confidence levels based on historical data. This method converts the uncertain PV energy from stochastic to deterministic to be included in the optimization model. The concept of demand response is extended and mathematically modeled using a linear function based on the quantized flexibility interval of multi-energy load demand. As a result, the overall model is formulated as a mixed-integer linear program, which can be effectively solved by the commercial solvers. The proposed model is implemented on two typical summer and winter days for various cases. Results of case studies show the important benefits for maximum PV utilization, energy efficiency and economic system operation. Moreover, the influence of the different confidence levels of PV power and effectiveness of IDR in the stochastic circumstances are addressed in the optimization-based operation.