Using the sub-seasonal to seasonal forecast model of Beijing Climate Center, several key physical parameters are perturbed by the Latin hypercube sampling method to find a better configuration for representation of Madden-Julian oscillation (MJO) in the free-run simulation. We find that although model simulation is especially sensitive to some parameters, there are overall no significant linear relationships between model skill and any one of the parameters, and the optimum performance can be obtained by combined perturbations of multiple parameters. By optimization, MJO's spectrum, intensity, spatial structure and propagation, as well as the mean state and variance, are all improved to some extent, suggesting the correspondence and interrelation of model's performances in simulating different characteristics of MJO. Further, several sets of initialized hindcasts using the optimized parameters are conducted, and their results are compared with the hindcasts using only improved initial conditions. We show that with an optimized model, the forecast of MJO beyond 3-week lead time is not improved, and the maximum useful skill is only slightly increased, implying that a decrease of model error does not always translate into an increase of forecast skill at all lead time. However, the skill is obviously enhanced during lead times of 2-3 weeks for forecasts in most seasons and initial phases except for a few cases. Particularly, the deficiency in forecasting MJO's propagation from the Indian Ocean to the Pacific is relieved, further highlighting the positive contribution of reducing model error compared to previous work that only reduced initial condition error. In this study, we also show benefits of multi-scheme ensemble strategy in describing uncertainties of model error and initial condition error and thus improving MJO forecast.