The Dual-Bubble-Size (DES) model is an extension of the Energy-Minimization Multi-Scale (EMMS) model, which was originally proposed for gas-solid fluidization to gas-liquid systems. This model is featured by a stability condition in addition to conservative equations for two bubble classes. The stability condition is mathematically formulated as a minimization tendency of microscale energy dissipation and physically reflects the compromise between different dominant mechanisms. This work attempts to use the Simulated Annealing (SA) method to solve the nonlinear optimization problem of the model. It is found that the SA method could greatly reduce the computational cost while capturing the jump change of gas holdup more accurately for the DBS model. The jump change reflects the transition from homogeneous and transitional regimes to heterogeneous regime, and essentially arises from the inflection of the trajectory of global minimum point within the space of the three structure parameters. The DBS model then is extended to Triple-Bubble-Size (TBS) and Multiple-Bubble-Size (MBS) models by introducing three or more bubble classes. We find that the TBS and MBS model prediction is reduced to that of DBS model and the two characteristic bubble classes are distinct in the calculation, even though the gas is resolved into three or more bubble classes. This implies that gas-liquid flow in bubble columns are essentially dominated by two bubble dasses rather than multiple bubble classes, and the DBS model may be an intrinsic model for reflecting the compromising mechanisms in the system and thereby describing the system evolution and gas-liquid interaction. Its integration with computational fluid dynamics (CFD) simulation may offer a more reasonable way to model the complex gas-liquid flow in bubble columns.