Civil structures and its energy dissipation devices often exhibit hysteretic and frequency-depended behavior.In this regard,the fractional derivative model can be used for frequency-dependent force-displacement relationship with very few parameters.The paper presents a multi harmonic balance method for determining steady-state response of a hysteretic system endowedwith fractional elements.Specifically,the Fourier series expansion both for the excitation and response and utilizing the multi harmonic balancing on both sides of the equation of motion,leads to a relationship between the Fourier coeficients of the excitationand that of the response.Next,recasting the hysteretic differential equation into a residue form and projecting it on the Fourier basisleads an algebraic equation with unknown response Fourier coeficients.The Levenberg-Marquardt algorithm is used to obtain theunknowns.Pertinent numerical examples including a softening and hardening Bouc-Wen system with different fractional order demonstrate the accuracy of the proposed method.