  We present a simulation-based inference (SBI) cosmological analysis of cosmic
shear two-point statistics from the fourth weak gravitational lensing data
release of the ESO Kilo-Degree Survey (KiDS-1000). KiDS-SBI efficiently
performs non-Limber projection of the matter power spectrum via 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, incorporating shear
measurement and redshift calibration uncertainties, as well as angular
anisotropies due to variations in depth and point-spread function. To enable
direct comparison with standard inference, we limit our analysis to
pseudo-angular power spectra. The SBI is based on sequential neural likelihood
estimation to infer the posterior distribution of spatially-flat $\Lambda$CDM
cosmological parameters from 18,000 realisations. We infer a mean marginal of
the growth of structure parameter $S_{8} \equiv \sigma_8 (\Omega_\mathrm{m} /
0.3)^{0.5} = 0.731\pm 0.033$ ($68 \%$). We present a measure of goodness-of-fit
for SBI and determine that the forward model fits the data well with a
probability-to-exceed of $0.42$. For fixed cosmology, the learnt likelihood is
approximately Gaussian, while constraints widen compared to a Gaussian
likelihood analysis due to cosmology dependence in the covariance. Neglecting
variable depth and anisotropies in the point spread function in the model can
cause $S_{8}$ to be overestimated by ${\sim}5\%$. Our results are in agreement
with previous analysis of KiDS-1000 and reinforce a $2.9 \sigma$ tension with
constraints from cosmic microwave background measurements. This work highlights
the importance of forward-modelling systematic effects in upcoming galaxy
surveys.
