Solving this non-holonomic system using Dirac-Bergmann theory I have read in some books and articles that the Dirac-Bergmann procedure to deal with constraints in phase space does not care about holonomic and Non-holonomic constraints, but I've been unable to find a single example. So, I wanted to test that assertion by solving a simple non-holonomic system.
Consider a particle ($m=1$) subject to the non-holonomic constraint $$\phi_{1}=\dot{y}-z\dot{x}=0.\tag{1}$$
The Lagrangian of the system is the standard one
$$L=\frac{1}{2}\left(\dot{x}^{2}+\dot{y}^{2}+\dot{z}^{2}\right)-V(\mathbf{r}),\tag{2}$$
and it is non-singular as the momentum can be found to be
$$p_{i}=\frac{\partial L}{\partial\dot{x}_{i}}=\dot{x}_{i}.$$
In phase-space, the dynamic is given by the Hamiltonian
$$H=\frac{1}{2}\left(p_{x}^{2}+p_{z}^{2}+p_{z}^{2}\right)+V(\mathbf{r})$$
constrained to obey $\phi_{1}=\dot{p_{y}}-z\dot{p_{z}}=0$. The time evolution is obtained using the Dirac bracket
$$\dot{F}=\left\{ F,H\right\} _{D}=\left\{ F,H\right\} -\sum_{i,j}\left\{ F,\phi_{i}\right\} \left(M_{ij}\right)^{-1}\left\{ \phi_{j},H\right\}, $$
where the Matrix of constraint has the following entries
$$M_{ij}=\left\{ \phi_{i},\phi_{j}\right\}.$$
Now, with only one constraint, the matrix only has one element, and since $\left\{ \phi_{1},\phi_{1}\right\} =0,$ the matrix is non-invertible and there is no Dirac Bracket.
I tried to remedy this in the usual way of the Dirac-Bergmann theory by introducing a second constraint
$$\phi_{2}=\left\{ \phi_{1},H\right\} \approx0$$
But the equations of motion that come from the Dirac bracket do not coincide with the ones from Lagrangian mechanics (analytical mechanics of discrete systems by Rosenberg, p 257).
So, given the above Hamiltonian and the constraint, how can the correct equation of motion be found?
 A: TL;DR: First of all, a constraint in the Dirac-Bergmann analysis of the Hamiltonian formulation cannot depend on velocities/dots. In contrast, OP's constraint (1) depends on velocities/dots.
Longer explanation: Well, there are at least 2 different Lagrangian formulations of non-holonomic constraints, cf. Ref. 1 and this Phys.SE post:

*

*One using a stationary action principle (SAP) but incompatible with d'Alembert's principle (DAP).


*One using Chetaev's equations, compatible with DAP but incompatible with SAP.
Since we need SAP to perform a (possible singular) Legendre transform to the Hamiltonian formulation, we will assume option 1.
Let us perform the first steps in the Dirac-Bergmann analysis of OP's example$^1$ (and leave it to the reader to complete it):
$$\begin{align} 
L~=~&L_0+\lambda\phi_1, \cr 
L_0 ~=~& \frac{1}{2}(\dot{x}^2+\dot{y}^2+\dot{z}^2)-V(x,y,z), \cr 
0~\approx~&\phi_1~=~\dot{y}-z\dot{x}. \end{align} \tag{A}$$
The Lagrangian momenta are
$$\begin{align}p_x~=~&\frac{\partial L}{\partial \dot{x}}~=~\dot{x}-\lambda z, \cr
p_y~=~&\dot{y}+\lambda, \cr
p_z~=~&\dot{z}, \cr
p_{\lambda}~=~&0. \end{align} \tag{B}$$
There is 1 primary constraint
$$0~\approx~\phi_2~=~p_{\lambda}.  \tag{C}$$
The original Hamiltonian becomes
$$\begin{align} H_0~=~&p_x\dot{x}+p_y\dot{y}+p_z\dot{z}+p_{\lambda}\dot{\lambda}-L\cr
~=~&\frac{1}{2}((p_x+\lambda z)^2+(p_y-\lambda)^2+p_z^2)+V.\end{align} \tag{D}$$
Note that $\lambda$ is not a Lagrange multiplier in the Hamiltonian formulation.
There is a secondary constraint
$$0~\approx~\phi_3~=~\{H_0,\phi_2\}~=~z(p_x+\lambda z)+\lambda-p_y. \tag{E} $$
There is a tertiary constraint
$$0~\approx~\phi_4~=~\{\phi_3,H_0\}~=~(p_x+2\lambda z)p_z -z\frac{\partial V}{\partial x}+\frac{\partial V}{\partial y}. \tag{F} $$
Here we stopped calculating, but it seems that the Hamiltonian formulation of OP's example becomes overconstrained.
References:

*

*P. Mann, Lagrangian and Hamiltonian dynamics, 2018; eqs. (8.2.4) & (8.2.4).


*R.M. Rosenberg, Analytical Dynamics of Discrete Systems, 1977; eq. (15.4.5).
--
$^1$ Ref. 2 uses option 2.
