How to get DAE representation from equations?

Update2: Thank you for all of your answers. I would like to clarify the question as I did a bit of research on DAE representations and their applications in software. There are 2 representations used in MATLAB:

The DAEs must be index 1 or 0 in both cases.

Both forms can be created using MATLAB Symbolic Toolbox. However, the consistent initial conditions function works only with the implicit form.

So my question is: how to get these DAE representations in Mathematica? I prefer to use Mathematica in symbolic computations, but in this case I may have to use MATLAB Symbolic.

Source of pictures: https://en.wikipedia.org/wiki/Differential-algebraic_system_of_equations

Update1: In the documentation I found methods to get the processed variables and equations. However, I still don't know how to use them to get the form that can be solved in MATLAB. The names of the variables are confusing and I'm not sure why there are more variables than equations.

I mentioned NonlinearStateSpaceModel because the output of this function can be used directly in MATLAB, but it is only for ODEs.

ClearAll["Global*"]

px := PX[t]
vx := D[PX[t], t]
ax := D[D[PX[t], t], t]
py := PY[t]
vy := D[PY[t], t]
ay := D[D[PY[t], t], t]

m := 1
g := -10
r := 1

Fx := z[t]*px
Fy := z[t]*py + m*g

vars := {PX, PY, z}

eqns = {r^2 == px^2 + py^2,
m*ax == Fx,
m*ay == Fy}

state = First@
NDSolveProcessEquations[{eqns, PX == 1, PY == 0}, vars, t,
Method -> {"IndexReduction" -> Automatic}]
residual = state["NumericalFunction"]["FunctionExpression"]
MatrixForm[residual[]]
MatrixForm[residual[]]

Original question: I have the following code which describes a simple pendulum using Newton's law of motion. This is a DAE (differential-algebraic equation) with higher order than one. NDSolve[] can solve the DAE with the settings.

Instead of numerical solving, I would like to get an order-reduced symbolic representation of the DAE: something similar like on this page: https://www.mathworks.com/help/matlab/math/solve-differential-algebraic-equations-daes.html

How to do it in Mathematica? NonlinearStateSpaceModel works really well with ODEs, but I cannot find such a nice function for DAEs.

ClearAll["Global*"]

px := PX[t]
vx := D[PX[t], t]
ax := D[D[PX[t], t], t]
py := PY[t]
vy := D[PY[t], t]
ay := D[D[PY[t], t], t]

m := 1
g := -10
r := 1

Fx := z[t]*px
Fy := z[t]*py + m*g

eqns = {r^2 == px^2 + py^2,
m*ax == Fx,
m*ay == Fy}

sol = NDSolve[{eqns, PX == 1, PY == 0, PX' == 0,
PY' == 0}, {PX, PY, z}, {t, 0, 5},
Method -> {"IndexReduction" -> Automatic}]

ParametricPlot[Evaluate[{PX[t], PY[t]} /. sol], {t, 0, 5}]
ParametricPlot[Evaluate[{PX[t], z[t]} /. sol], {t, 0, 5}]
• Have you looked at the section Index Reduction for DAEs in this documentation article? Jan 7 '20 at 2:09
• The plots are irrelevant to your question, no? Jan 7 '20 at 4:04
• It seems the function expression residual is the $f$ in the implicit form $f=0$ used by the IDA method in NDSolve, like MATLAB's ode15i (which may or may not use IDA under the hood). I still can't answer the question about why there are more variables than equations. Jan 8 '20 at 16:58

To get the first-order state-space equations you can use an internal function that NonlinearStateSpaceModel uses under the hood. NonlinearStateSpaceModel cannot handle this at the moment because it has algebraic equations, but the internal function has no problem with such equations.

{{f, h, e}, x} = ControlDEqns`nonaffinestatespaceForm[
eqns, {PX[t], PY[t], z[t]}, {}, {PX[t], PY[t], z[t]}, t, #[[1 ;; 2]] &,
DescriptorStateSpace -> True];

The above will give equations of the form $$e.x'==f$$.

$$h$$ is the outputs which I specified as $$\{PX[t], PY[t], z[t]\}$$.  