ModelingToolkit.jl ODEProblem throwing error when differential term not isolated to the LHS

I am constructing an ODESystem by connecting 4 different systems.
I recently changed one of the systems, I want it to be a 2-DOF planar arm with equations

M(\theta)\ddot{\theta} + C(\theta,\dot{\theta})\dot{\theta} = u

M(\theta) = \begin{pmatrix} I_1+I_2+m_2l_1^2+2m_2l_1\bar{l}_2\cos(\theta_2) & I_2+m_2l_1\bar{l}_2\cos(\theta_2) \\ I_2+m_2l_1l_2\cos(\theta_2) & I_2 \end{pmatrix}

C(\theta,\dot{\theta}) = m_2l_1l_2\sin(\theta_2)\begin{pmatrix}-2\dot{\theta}_2 & -\dot{\theta}_2\\ \dot{\theta}_1 & 0 \end{pmatrix}+ \begin{pmatrix} D_1 & 0 \\ 0 & D_2 \end{pmatrix}

I create this system with the following code

num_inputs=2
@parameters t, m1, m2, l1, l2, l2b, I1, I2, D1, D2
dt = Differential(t)

@variables x[1:num_inputs](t)
@variables u[1:num_inputs](t)
@variables o(t)

# equations for the plant
eqs = vcat(
      M(x[2],I1,I2,m2,l1,l2,l2b)*dt.(dt.(x)) .+ Cf(x[1],x[2],dt(x[1]),dt(x[2]),m2,l1,l2,D1,D2)*dt.(x) .~ u,
      o .~ x[1]
  )    
return ODESystem(eqs, t; name=name,defaults = Dict([m1,m2,l1,l2,l2b,I1,I2,D1,D2] .=> defaultP))  

The system generation works well

Model NLP with 3 equations
States (5):
  x[1](t)
  x[2](t)
  u[1](t)
  u[2](t)
  o(t)
Parameters (8):
  D1 [defaults to 0.05]
  l2 [defaults to 0.35]
  l1 [defaults to 0.3]
  l2b [defaults to 0.21]
  m2 [defaults to 2.5]

the equations of the system are

equations(psys)
3-element Vector{Equation}:
 (I2 + l1*l2b*m2*cos(x[2](t)))*Differential(t)(Differential(t)(x[2](t))) + (D1 - 2l1*l2*m2*Differential(t)(x[2](t))*sin(x[2](t)))*Differential(t)(x[1](t)) + (I1 + I2 + m2*(l1^2) + 2l1*l2b*m2*cos(x[2](t)))*Differential(t)(Differential(t)(x[1](t))) - l1*l2*m2*(Differential(t)(x[2](t))^2)*sin(x[2](t)) ~ u[1](t)
 D2*Differential(t)(x[2](t)) + I2*Differential(t)(Differential(t)(x[2](t))) + (I2 + l1*l2*m2*cos(x[2](t)))*Differential(t)(Differential(t)(x[1](t))) + l1*l2*m2*(Differential(t)(x[1](t))^2)*sin(x[2](t)) ~ u[2](t)
 o(t) ~ x[1](t)

Then I connect this system to the other 3 systems, this works without issues. Then I structural_simplify the system which works again

sys = structural_simplify(all)
Model FFFBsys with 8 equations
NLP₊x[2](t)
  NLP₊x[1](t)
  intError(t)
  eI(t)
  NLP₊xˍt[2](t)
⋮
Parameters (268):
  NLP₊D1 [defaults to 0.05]
  NLP₊l2 [defaults to 0.35]
  NLP₊l1 [defaults to 0.3]
  NLP₊l2b [defaults to 0.21]
  NLP₊m2 [defaults to 2.5]
⋮
Incidence matrix:8×16 SparseArrays.SparseMatrixCSC{Num, Int64} with 26 stored entries:
⠑⢄⠀⠀⢠⠀⠑⠀
⠛⠀⠛⢄⢫⠀⠀⢉

the simplified system has the following equations:
equations(sys)

8-element Vector{Equation}:
 Differential(t)(NLP₊x[2](t)) ~ NLP₊xˍt[2](t)
 Differential(t)(NLP₊x[1](t)) ~ NLP₊xˍt[1](t)
 Differential(t)(intError(t)) ~ 0.5(e(t)^2)
 Differential(t)(eI(t)) ~ e(t)
 0 ~ NLP₊l1*NLP₊l2*NLP₊m2*(Differential(t)(NLP₊x[2](t))^2)*sin(NLP₊x[2](t)) + NLP₊u[1](t) - (NLP₊I2 + NLP₊l1*NLP₊l2b*NLP₊m2*cos(NLP₊x[2](t)))*Differential(t)(NLP₊xˍt[2](t)) - (NLP₊D1 - 2NLP₊l1*NLP₊l2*NLP₊m2*Differential(t)(NLP₊x[2](t))*sin(NLP₊x[2](t)))*Differential(t)(NLP₊x[1](t)) - (NLP₊I1 + NLP₊I2 + NLP₊m2*(NLP₊l1^2) + 2NLP₊l1*NLP₊l2b*NLP₊m2*cos(NLP₊x[2](t)))*Differential(t)(NLP₊xˍt[1](t))
 0 ~ MLP₊o[2](t) - NLP₊D2*Differential(t)(NLP₊x[2](t)) - (NLP₊I2 + NLP₊l1*NLP₊l2*NLP₊m2*cos(NLP₊x[2](t)))*Differential(t)(NLP₊xˍt[1](t)) - NLP₊I2*Differential(t)(NLP₊xˍt[2](t)) - NLP₊l1*NLP₊l2*NLP₊m2*(Differential(t)(NLP₊x[1](t))^2)*sin(NLP₊x[2](t))
 Differential(t)(pid₊int[1](t)) ~ 0.1pid₊u[1](t)
 Differential(t)(pid₊mv[1](t)) ~ pid₊u[1](t) - 0.1pid₊mv[1](t)

The problem arises when I try to define and ODEProblem from the structurally simplified system, I get the following error:

InvalidSystemException: The derivative variable must be isolated to the left-hand side of the equation like `Differential(t)(NLP₊x[4](t)) ~ ...`. You may want to use `structural_simplify` or the DAE form.
Got 0 ~ NLP₊l1*NLP₊l2*NLP₊m2*(NLP₊x[4](t)^2)*sin(NLP₊x[2](t)) + NLP₊u[1](t) - (NLP₊I2 + NLP₊l1*NLP₊l2b*NLP₊m2*cos(NLP₊x[2](t)))*Differential(t)(NLP₊x[4](t)) - (NLP₊D1 - 2NLP₊l1*NLP₊l2*NLP₊m2*NLP₊x[4](t)*sin(NLP₊x[2](t)))*NLP₊x[3](t) - (NLP₊I1 + NLP₊I2 + NLP₊m2*(NLP₊l1^2) + 2NLP₊l1*NLP₊l2b*NLP₊m2*cos(NLP₊x[2](t)))*Differential(t)(NLP₊x[3](t)).

I guess the more general issue is: is there a way of defining an ODESystem when the equations don’t have an isolated differential on the LHS?
I have seen some similar issues that lead to this issue Automatically move differentials to the left side by ChrisRackauckas · Pull Request #1164 · SciML/ModelingToolkit.jl · GitHub
Is there a solution now implemented? Is there a workaround otherwise?

Thanks for the help!

Looks like your example code is missing the definition of M, so I did not test this, but you may need to use ODAEProblem instead of ODEProblem after simplifying.

https://mtk.sciml.ai/stable/mtkitize_tutorials/modelingtoolkitize_index_reduction/