The aim of this paper is to prove the existence of optimal shapes in bilinear parabolic optimal control problems. We consider a parabolic equation \(\partial _tu_m-\Delta u_m=f(t,x,u_m)+mu_m\) . The set of admissible controls is given by \(A=\{m\in L^\infty \,, m_-\leqq m\leqq m_+{\text { almost everywhere, }}\int _\Omega m(t,\cdot )=V_1(t)\}\) , where \(m_\pm =m_\pm (t,x)\) are two reference functions in \(L^\infty ({(0,T)\times {\Omega }})\) , and where \(V_1=V_1(t)\) is a reference integral constraint. The functional to optimise is \(J:m\mapsto \iint _{(0,T)\times {\Omega }} j_1(u_m)+\int _{\Omega }j_2(u_m(T))\) . Roughly speaking, we prove that, if \(j_1\) and \(j_2\) are non-decreasing and if one is increasing, then any solution of \(\max _A J\) is bang-bang: any optimal \(m^*\) writes \(m^*=\mathbb {1}_E m_-+\mathbb {1}_{E^c}m_+\) for some \(E\subset {(0,T)\times {\Omega }}\) . From the point of view of shape optimization, this is a parabolic analog of the Buttazzo-Dal Maso theorem in shape optimisation. The proof is based on second-order criteria and on an approximation-localisation procedure for admissible perturbations. This last part uses the theory of parabolic equations with measure data.