Solving systems of Boolean equations is a fundamental task in symbolic computation and algebraic cryptanalysis, with wide-ranging applications in cryptography, coding theory, and formal verification. Among existing approaches, the Boolean Characteristic Set (BCS) method[1] has emerged as one of the most efficient algorithms for tackling such problems. However, its performance is highly sensitive to the ordering of variables, with solving times varying drastically under different orderings for fixed variable counts n and equations size m. To address this challenge, this paper introduces a novel optimization framework that synergistically integrates machine learning (ML)-based time prediction with simulated annealing (SA) to efficiently identify high-performance variables orderings. Weconstruct a dataset comprising variable frequency spectrum X and corresponding BCS solving time t for benchmark systems(e.g., n = m = 28). Utilizing this data, we train an accurate ML predictor ft(X) to estimate solving time for any given variables ordering. For each target system, ft serves as the cost function within an SA algorithm, enabling rapid discovery of low-latency orderings that significantly expedite subsequent BCS execution. Extensive experiments demonstrate that our method substantially outperforms the standard BCS algorithm[1], Gröbner basis method [2] and SAT solver[3], particularly for larger-scale systems(e.g., n = 32). Furthermore, we derive probabilistic time complexity bounds for the overall algorithm using stochastic process theory, establishing a quantitative relationship between predictor accuracy and expected solving complexity. This work provides both a practical acceleration tool for algebraic cryptanalysis and a theoretical foundation for ML-enhanced combinatorial optimization in symbolic computation.