This paper presents an approximate, frequency-domain approach to modelling complex structures (CSs) with localised nonlinearities, designated here as complex sub-systems. This consistent approach to modelling CSs presented here aims to improve the computational efficiency, which in cases when nonlinearities are included, is a major problem when CSs with many degrees of freedom are modelled. This approach proposes sub-structuring of the CS into its linear and nonlinear parts in the first stage. Classical reduction of the linear part and the nonlinear model reduction and the polynomial approximation for the nonlinear part are employed in the second stage to decrease the overall number of degrees of freedom. Finally, an additional, well-suited harmonic-balance and describing-function-based approximation is used for the nonlinear part, introducing the multi-coordinate describing functions (MCDFs) and the multi-coordinate describing-function matrix (MCDFM). Together with the matrices of the linear part of the localised nonlinearity, the MCDFM forms the so-called harmonic nonlinear super-model (HNSM). The HNSM introduced is well-suited for use with FRF coupling in the frequency domain. Two numerical case studies as well as an experimental case study showed that this approach is suitable for the steady-state vibration of CSs with localised nonlinearities, while at the same time, the efficient approach makes it possible to perform parametric analyses. It is shown that, with some restrictions, an optional reconstruction is also possible, which makes this approach even more efficient.