An algorithm for imaginary time evolution of a fermionic projected entangled pair state with ancillas from infinite temperature down to a finite temperature state is presented. As a benchmark application, it is applied to spinless fermions hopping on a square lattice subject to p-wave pairing interactions. With a tiny bias it allows one to evolve the system across a high-temperature continuous symmetry-breaking phase transition.