High-order stochastic simulation overcomes the limitation of conventional geostatistical simulation methods in reproducing complex spatial patterns in natural phenomena. It also has the data-driven characteristic, an advantage over the multiple-point simulations that are commonly training-image driven. High-order stochastic simulation, in contrast to assuming a Gaussian distribution as in conventional geostatistical simulation methods, the probability distribution of the underlying random field is approximated by an orthogonal polynomial expansion series. The approximation takes place as the high-order spatial statistics of the available data, or a training image, are utilized to derive the coefficients of the polynomial expansion series. Under the conceptual framework of high-order geostatistical simulation, high-order spatial statistics characterize statistical interactions among multiple points and act as a mathematical tool capturing the complex spatial patterns of natural attributes, similar in a way as variograms to capture the two-point statistics in conventional simulation methods. As a result, high-order simulation has the capacity to reproduce complex spatial patterns from the available data. On the other hand, a challenge is to address its numerical stability, in particular, for approximating a conditional probability distribution along the simulation procedure. This challenge becomes prominent when statistical conflicts exist between the actual samples and a training image, leading to negative values in probability related to numerical instability. A mathematical solution is proposed to address the above-mentioned challenge in the context of mathematical optimization. Specifically, a new approximation by polynomial series is achieved by approaching the high-order statistics of the available data while simultaneously stabilizing the solution via semidefinite programming. This approach guarantees the positiveness of the approximation of a distribution and eliminates the possible numerical instability in high-order simulations. Examples demonstrate the efficacy of the proposed method.