This work proposes an iterative identification algorithm for polynomial NARX models with several inputs. An optimal input design procedure is used to obtain persistent excitation for all regressors taking into account actuator and safety constraints. Through this procedure and the iterative increase of the polynomial degree, systems of unknown nonlinear structure can be approximated. The paper presents the procedure and results using a numerical and experimental example – the air-path of a combustion engine.