We perform for the first time full simulation-based inference on the Lyman-$\alpha$ forest 1D power spectrum. In particular, we consider the prediction of the Lyman-$\alpha$ forest $P_{\rm 1D}(k)$ at $2.0<z<3.5$ from the cosmological hydrodynamic simulations run with the and galaxy formation models. We train a normalizing flow to perform neural posterior estimation of two cosmological parameters (${\Omega }_m$ and ${\sigma }_8$) and four astrophysical parameters parametrizing supernova and AGN feedback. When training and testing the neural network on the same baryon physics model, the posterior distributions of the cosmological parameters are found to be in excellent agreement with the true parameters values (within $10\%$ deviations in $\gtrsim 75\%$ and $\gtrsim 90\%$ of the cases for ${\Omega }_m$ and ${\sigma }_8$, and a precision better than $10\%$ in both), while the astrophysical parameters converge to the prior mean due to the limited probed volume. When training on one model and testing on the other (e.g., training on and testing on , or viceversa), the performance is significantly worse, both in accuracy and in precision, resulting in a $\sim 10\%$ positive bias on the predicted values for ${\sigma }_8$. We show that a multi-domain training based on the combination of simulations from both models recovers unbiased constraints, offering an effective solution to cope with the complex problem of the lack of convergence in the predictions from different galaxy formation models. This study represents a promising way forward to constrain cosmology and fundamental physics with the Lyman-$\alpha$ forest with artificial intelligence.