The Gibbs-ensemble Monte Carlo method and the Gibbs-Duhem integration scheme are adapted for the simulation of the phase equilibrium of a one-component fluid confined in random porous media. The validity of these methods in the case of rigid porous samples is established by comparing our results with those obtained previously from a series of adsorption isotherms. It is shown that the Gibbs-ensemble and Gibbs-Duhem integration methods significantly improve the efficiency of the simulation of these systems. Such a gain in efficiency allowed us to carry out a systematic investigation of the influence of several characteristics of disordered porous solids (e.g., porosity, pore size distribution, and solid-fluid interaction) on the phase behavior of the confined fluid. Rich phase behaviors have been observed, e.g., multiple fluid-fluid phase transitions and an extreme sensitivity of phase diagram on the microscopic structure of the porous samples. Efforts were devoted to understanding the origins of such rich behavior by analyzing the simulation results in considerable detail.