In the derivation of the generating function of the Gaudin Hamiltonians with boundary terms we follow the same approach used previously in the rational case, which in turn was based on Sklyanin's method in the periodic case. Our derivation is centered on the quasiclassical expansion of the linear combination of the transfer matrix of the XXZ Heisenberg spin chain and the central element, the so-called Sklyanin determinant. The corresponding Gaudin Hamiltonians with boundary terms are obtained as the residues of the generating function. By defining the appropriate Bethe vectors which yield strikingly simple off-shell action of the generating function, we fully implement the algebraic Bethe ansatz, obtaining the spectrum of the generating function and the corresponding Bethe equations. * The so-called rational sℓ(2) Gaudin model was introduced by Gaudin himself as a model of "long-range" interacting spins in a chain [1]. This model was extended to any simple Lie algebra, with arbitrary irreducible representation at each site of the chain [2,3]. Sklyanin studied the rational model in the framework of the quantum inverse scattering method [4,5] using the sℓ(2) invariant classical r-matrix [6]. A generalization of these results to all cases when skewsymmetric r-matrix satisfies the classical Yang-Baxter equation [7] was relatively straightforward [8,9]. Therefore, considerable attention has been devoted to Gaudin models corresponding to the the classical r-matrices of simple Lie algebras [10][11][12][13][14][15] and Lie superalgebras [16][17][18][19][20].Hikami, Kulish and Wadati showed that the quasi-classical expansion of the transfer matrix of the periodic chain, calculated at the special values of the spectral parameter, yields the Gaudin Hamiltonians [21,22]. Hikami showed how the quasi-classical expansion of the transfer matrix, calculated at the special values of the spectral parameter, yields the Gaudin Hamiltonians in the case of non-periodic boundary conditions [23]. Then the ABA was applied to open Gaudin model in the context of the the Vertex-IRF correspondence [24][25][26]. Also, results were obtained for the open Gaudin models based on Lie superalgebras [27]. An approach to study the open Gaudin models based on the classical reflection equation [28][29][30] and the non-unitary r-matrices was developed recently, see [31,32] and the references therein. For a review of the open Gaudin model see [33]. Progress in applying Bethe ansatz to the Heisenberg spin chain with non-periodic boundary conditions compatible with the integrability of the quantum systems [34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51] had recent impact on the study of the corresponding Gaudin model [52,53]. The so-called T − Q approach to implementation of Bethe ansatz [39,40] was used to obtain the eigenvalues of the associated Gaudin Hamiltonians and the corresponding Bethe ansatz equations [54]. In [51] the off-shell action of the generating function of the trigonometric Gaudin Hamiltonians with boundary terms o...