In recent years, Deep Neural Networks (DNNs) have emerged as a new framework for solving Partial Differential Equations (PDEs). The most widely used approach, Physics-Informed Neural Networks (PINNs), minimizes the L2 norm of the residual. This formulation imposes restrictive regularity requirements on the solution. To address these limitations, several variational formulations have been proposed in the literature. The Deep Ritz method minimizes an energy functional associated to the problem and Variational PINNs minimize a weak residual projected onto a finite-dimensional test space. However, these functionals fail to satisfy an interpolation property, which results in training instabilities due to stochastic integration errors. The First-Order System Least Squares (FOSLS) methodology recasts higher-order PDEs as a first-order system for which the L2 norm of the residual is variationally consistent. Unlike the aforemetioned methods, the FOSLS functional provides inherent variance reduction that stabilizes training. However, standard FOSLS functionals yield residual–error norm bounds with parameter-dependent constants. For transmission problems governed by the Poisson equation, we introduce a FOSLS functional that is robust with respect to a physically consistent energy norm of the error. A key contribution is our novel definition of the energy norm. Numerical experiments and theoretical results show that the proposed method is significantly less sensitive to stochastic quadrature errors than the Deep Ritz method and remains effective under coarse quadrature rules. Problems with discontinuous coefficients illustrate the advantages of the robust FOSLS functional and the adaptive basis approach with respect to existing methods such as PINNs, VPINNs, and Deep Ritz.