Mistake with your approach
Although your transformations should preserve the tangency conditions, this does not mean that the points of tangency are fixed points. In fact they can not be. Think about it this way: if you have two points and the tangent directions in these points, then these uniquely define the point where the tangents intersect. If all of this is to remain fixed, you have three fixed points and the only affine transformation with three non-collinear fixed points is the identity. Ergo your points of tangency can't be fixed.
Family of conics
So what's your family of conics? You already know one such conic, namely the unit circle:
$$x^2+y^2-1=0$$
Another one is the degenerate ellipse of zero minor semiaxis which passes through the two points of tangency:
$$(x+y-1)^2=x^2+y^2-2x-2y+2xy+1=0$$
The family of conics (or pencil of conics) would be the set of all linear combinations of these two. Introducing a parameter $p$ you could write them as
$$(x^2+y^2-1)+2p(x+y-xy-1)=0$$
or something like that.
Describing desired transformations
In order to find your affine transformation, start by writing this as a matrix:
$$(x,y,1)\cdot\begin{pmatrix}
1 & -p & p \\
-p & 1 & p \\
p & p & -1-2p
\end{pmatrix}\cdot\begin{pmatrix}
x\\y\\1\end{pmatrix}=0$$
Applying an affine transformation to a conic means conjugating the matrix of the conic with the matrix of the inverse transformation. If we want to apply the forward transformation instead, we can say that this forward transformation should turn the conic for parameter $p$ back into the unit circle:
$$\begin{pmatrix}
A & D & 0 \\
B & E & 0 \\
C & F & 1
\end{pmatrix}\cdot\begin{pmatrix}
1 & -p & p \\
-p & 1 & p \\
p & p & -1-2p
\end{pmatrix}\cdot\begin{pmatrix}
A & B & C \\
D & E & F \\
0 & 0 & 1
\end{pmatrix}=q\cdot\begin{pmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & -1
\end{pmatrix}$$
The variable $q$ on the right hand side is required because any multiple of that matrix will describe the same conic (so the matrix is homogeneous). Any non-zero value for $q$ is acceptable.
The space of all solutions
This matrix equation gives you nine scalar equations, but due to the matrix being symmetric, only six of them are actually distinct. And even among those there are interdependencies, so that for a given $p$ the set of solutions forms a two-dimensional variety. Roughly speaking you have two real degrees of freedom still to choose your transformation.
Specific solutions
To obtain a unique solution, you may want to require $A=E,B=D,C=F$, i.e. make your transformation fully symmetric with respect to swapping the $x$ and $y$ coordinates. With the above approach you end up with four quadratic equations for four variables $A,B,C,q$:
\begin{align*}
A^2p + B^2p - 2AB &= 0 \\
A^2 + B^2 - 2ABp &= q \\
2C^2p - 2C^2 - 4Cp + 2p + 1 &= q \\
(A+B)(C - Cp + p) &= 0 \\
\end{align*}
But depending on your application, other transformations might be of interest as well, e.g. this one here (for $p=-\frac12$ and with additional conditions $2A=E,C=F$):
$$\begin{pmatrix}
A & B & C \\
D & E & F \\
0 & 0 & 1
\end{pmatrix}=\begin{pmatrix}
\sqrt{\frac{5-2\sqrt3}{13}} & \sqrt{\frac{47+2\sqrt3}{13}} & -1 \\
2\sqrt{\frac{8+2\sqrt3}{13}} & 2\sqrt{\frac{5-2\sqrt3}{13}} & -1 \\
0 & 0 & 1
\end{pmatrix}$$
Non-ellipse solutions
Note that you can't expect to always find such a matrix for any $p$ you may choose. For example, for $p=1$ you get a parabola. So your family of conics contains parabolas and hyperbolas. You can't get these from affine transformations of a circle.
In my first versions of this post I had forgotten about the parameter $q$ above, leading me to conclude that no such transformation could exist. In order to offer some alternatives, I then discussed ways to use projective transformations instead. If you need those other conics as well, please look at the revision history of this post for details on that aspect.
Comparison to other answers
I also tried out Yves' answer for $p=-1$ in his formulas which corresponds to $p=\frac13$ in mine. I came up with $a=-\sqrt{1/6}, b=-\sqrt{1/3}$ and
$$\begin{pmatrix}
A & B & C \\
D & E & F \\
0 & 0 & 1
\end{pmatrix}=\begin{pmatrix}
-\sqrt{\frac32} & \frac{\sqrt3}2 & -\frac12 \\
-\sqrt{\frac32} & -\frac{\sqrt3}2 & -\frac12 \\
0 & 0 & 1
\end{pmatrix}$$
which looks more like $A=D,B=-E,C=F$ as additional constraint. I'm a bit surprised by this at the moment, but I might have made a mistake somewhere.