How to generate random matrices that satisfy Lyapunov equations

Question

I wish to generate a large number of matrix quadruples $(A,B,C,D)$ where $A\in \mathbb{R}^{n\times n}$, $B \in \mathbb{R}^{n \times 2}$, $C \in \mathbb{R}^{2 \times n}$ and $ D \in \mathbb{R}^{2 \times 2}$ that satisfy $$ A+A^T = -BB^T = -C^TC, \quad D^TD = I. $$ The condition is similar to Lyapunov equation with the unknown set to identity.

Any insight would be greatly appreciated.