Abstract：Based on the immune cell subsets division and immune life theory, an immune based multiobjective gene expression programming(GEP) algorithm for identifying the polynomial nonlinear auto regressive moving average with exogenous inputs(NARMAX) model is proposed. In the algorithm, the cell initialization method and gene operators are newly defined. Based on the multi-optimal-solutions characteristic of multiobjective optimization problems, the parameter performance of each subset is estimated by the optimal solutions of the subset in the total optimal solutions, and then the search direction for the next generation is determined. The simulation results show that the proposed algorithm performs well in the polynomial NARMAX modeling, and the structure and the parameters of the nonlinear system can be identified at the same time.