There are much research effort in the literature using Monte Carlo simulation (MCS) which is a direct and simple numerical method, however, it can be computationally very expensive as the governing dynamic equations of the system to be simulated for each random sample using the MCS. In this paper, polynomial meta-models based on the evolved group method of data handling (GMDH) neural networks are obtained to simply calculate the probability of failure in the MCS, instead of direct solution of dynamic equation of system. In this way, some input–output data consisting of uncertain parameters of system and controller parameters as inputs and probability of failure of some cost functions as output are used for training such GMDH-type neural networks which replace the very time consuming direct solution of dynamical systems during the MCS. A multi-objective genetic algorithm is also used for Pareto optimal design of PI and PID controllers for both first and second order uncertain system with time delays using methodology of this paper. The objective functions that are considered for such Pareto multi-objective optimization are namely, probability of failure of settling time (PTS) and probability of failure of overshoot (POS). The comparisons of the obtained results using the method of this paper with those obtained using direct method shows a significant reduction in computational time, whilst the accuracy is maintained.