Abstract: Components or systems fails are usually caused by multiple failure models in practical engineering. Copula functions are used to model complicated nonlinear corrections between variables and failure modes. To avoid error from subjective hypothesis at utmost, maximum entropy approach is employed for approximating probability density functions of extreme value responses, and the lower and upper bounds of probability of failure can be determined according to the approximated probability density functions. Finally, both simulation and maximum likelihood estimation approaches are used for determining the lower and upper bounds of parameters in selected Copula functions. An optimization model is presented to calculate the lower and upper bounds of system probability of failure. Numerical example has demonstrated that the proposed method is available for structural systems with multiple failure modes under mixed uncertainties.