Let us consider the system depicted in picture.
The position of the three points of interest can be easily found as:
$p_a = \begin{pmatrix} \epsilon \\ 0\end{pmatrix}$
$p_b = p_a + R_\alpha \begin{pmatrix} l_1 \\ 0\end{pmatrix}, \qquad R_\alpha = \begin{pmatrix} \cos \alpha & -\sin \alpha \\ \sin \alpha & \cos \alpha \end{pmatrix}$
$p_c = p_b + R_\alpha R_\beta \begin{pmatrix} l_2 \\ 0\end{pmatrix}, \qquad R_\beta = \begin{pmatrix} \cos \beta & -\sin \beta \\ \sin \beta & \cos \beta \end{pmatrix}$
Considering the joint variables:
$q = \begin{pmatrix} \epsilon, & \alpha, & \beta \end{pmatrix}^\top$
we can formulate the control problem of finding the control law to reduce asymptotically the error:
$\tilde{p} = p_d – p_c(q)$
between a fixed position $p_d \in \mathbb{R}^2$ such that $\dot{p}_d = 0$ and the result of the direct kinematics represented by the function $p_c(q): \mathbb{R}^3 \rightarrow \mathbb{R}^2$.
Let us consider a Lyapunov function:
$V(\tilde{p}) = \frac{1}{2} \tilde{p}^\top \tilde{p} \succ 0 \quad \forall \quad \tilde{p} \neq 0$
computing the time derivative we obtain:
$\dot{V}(\tilde{p}) = \tilde{p}^\top (-J(q)\dot{q}), \qquad J(q) =\frac{\partial p_c(q)}{\partial q}$
where $J(q)$ is the analytical Jacobian.
As the choice of the control action should enforce negative definiteness for the Lyapunov function time derivative, a possible choice can be immediately identified by: $\dot{q} = K J^\top(q) \left( p_d – p_c(q) \right)$
It is well known that kinematics singularities might affect the control results as the matrix $J(q) K J^\top(q)$ is responsible of the quadratic form: $\dot{V}(\tilde{p}) = -\tilde{p}^\top J(q) K J^\top(q) \tilde{p} \prec 0$
A good indicator of the singularity proximity is given by the manipulability computed as: $\mu =det(J(q)J^\top(q))$

Let us now convert our math results into runnable code with Sym2Code