KiDS-SBI is a simulation-based inference method used to analyze cosmic shear two-point statistics from the fourth weak gravitational lensing data release of the ESO Kilo-Degree Survey (KiDS-1000). The method efficiently performs non-Limber projection of the matter power spectrum using Levin's method and constructs log-normal random matter fields on the curved sky for arbitrary cosmologies, including effective prescriptions for intrinsic alignments and baryonic feedback. The forward model samples realistic galaxy positions and shapes based on the observational characteristics of KiDS-1000, incorporating shear measurement and redshift calibration uncertainties, as well as angular anisotropies due to variable survey depth and point-spread function variations. The analysis is limited to pseudo-angular power spectra as summary statistics. The SBI is based on neural density estimation of the likelihood with active learning to infer the posterior distribution of spatially-flat Λ CDM cosmological parameters from 18,000 realisations. The analysis infers a mean marginal of the growth of structure parameter S₈ ≡ σ₈(Ωₘ/0.3)^0.5 = 0.731 ± 0.033 (68%). The results agree with previous analysis of KiDS-1000 and reinforce a 2.9σ tension with early-Universe constraints from cosmic microwave background measurements. This work highlights the importance of forward-modelling systematic effects in upcoming galaxy surveys, such as Euclid, Rubin and Roman.KiDS-SBI is a simulation-based inference method used to analyze cosmic shear two-point statistics from the fourth weak gravitational lensing data release of the ESO Kilo-Degree Survey (KiDS-1000). The method efficiently performs non-Limber projection of the matter power spectrum using Levin's method and constructs log-normal random matter fields on the curved sky for arbitrary cosmologies, including effective prescriptions for intrinsic alignments and baryonic feedback. The forward model samples realistic galaxy positions and shapes based on the observational characteristics of KiDS-1000, incorporating shear measurement and redshift calibration uncertainties, as well as angular anisotropies due to variable survey depth and point-spread function variations. The analysis is limited to pseudo-angular power spectra as summary statistics. The SBI is based on neural density estimation of the likelihood with active learning to infer the posterior distribution of spatially-flat Λ CDM cosmological parameters from 18,000 realisations. The analysis infers a mean marginal of the growth of structure parameter S₈ ≡ σ₈(Ωₘ/0.3)^0.5 = 0.731 ± 0.033 (68%). The results agree with previous analysis of KiDS-1000 and reinforce a 2.9σ tension with early-Universe constraints from cosmic microwave background measurements. This work highlights the importance of forward-modelling systematic effects in upcoming galaxy surveys, such as Euclid, Rubin and Roman.