The Fuzzy Broad Learning System (Fuzzy BLS) is established by replacing the feature nodes of a Broad Learning System with the Takagi–Sugeno–Kang (TSK) fuzzy sub-systems. K-means algorithm is employed to cluster the input data so as to reduce computation complexity. And the parameters of a Fuzzy BLS are computed analytically by pseudoinverse. We investigate the learning algorithms of Fuzzy BLS comprehensively and apply them to nonlinear system identification in this paper: First of all, we develop an iterative learning algorithm for updating the weights in top layer and the weights connecting the fuzzy sub-systems to the enhancement nodes by gradient descent. Secondly, we analyze and prove the Lyapunov stability of Fuzzy BLS with this iterative algorithm. Then, we consider the fuzzy c-means for clustering input data in the part of fuzzy sub-systems, as well as randomly generated centers for Gaussian membership functions. There are several different learning algorithms due to the choice of clustering methods and calculating parameters by pseudoinverse or gradient descent iteratively, which are compared with each other in detail by system identification problems. It is concluded that the learning algorithms which calculate weights by pseudoinverse always outperform the ones that update them iteratively, no matter which clustering method is chosen. The fuzzy c-means, c-means and random centers each has its own merits in our experiments. In addition, Fuzzy BLS trained by the proposed algorithms demonstrates its superiority over the state-of-the-art neuro-fuzzy models in identifying nonlinear systems.