In methods like geminal-based approaches or coupled cluster that are solved using the projected Schrödinger equation, direct computation of the 2-electron reduced density matrix (2-RDM) is impractical and one falls back to a 2-RDM based on response theory. However, the 2-RDMs from response theory are not $N$-representable. That is, the response 2-RDM does not correspond to an actual physical $N$-electron wave function. We present a new algorithm for making these non-$N$-representable 2-RDMs approximately $N$-representable, i.e., it has the right symmetry and normalization and it fulfills the $P$-, $Q$-, and $G$-conditions. Next to an algorithm which can be applied to any 2-RDM, we have also developed a 2-RDM optimization procedure specifically for seniority-zero 2-RDMs. We aim to find the 2-RDM with the right properties which is the closest (in the sense of the Frobenius norm) to the non-$N$-representable 2-RDM by minimizing the square norm of the difference between this initial response 2-RDM and the targeted 2-RDM under the constraint that the trace is normalized and the 2-RDM, $Q$-matrix, and $G$-matrix are positive semidefinite, i.e., their eigenvalues are non-negative. Our method is suitable for fixing non-$N$-representable 2-RDMs which are close to being $N$-representable. Through the $N$-representability optimization algorithm we add a small correction to the initial 2-RDM such that it fulfills the most important $N$-representability conditions.