ModelingToolkitStandardLibrary.jl icon indicating copy to clipboard operation
ModelingToolkitStandardLibrary.jl copied to clipboard

Start `PlanarMechanics` module

Open AbdulrhmnGhanem opened this issue 2 years ago • 14 comments

Inspired by https://github.com/dzimmer/PlanarMechanics

I'll add all the basic components from the PlanarMechanics Modelica library. Are you OK with adding it to the standard library?

AbdulrhmnGhanem avatar Sep 23 '23 05:09 AbdulrhmnGhanem

Codecov Report

Attention: Patch coverage is 0.29412% with 339 lines in your changes are missing coverage. Please review.

Project coverage is 12.98%. Comparing base (35a41f7) to head (8838212). Report is 5 commits behind head on main.

:exclamation: Current head 8838212 differs from pull request most recent head 611a124. Consider uploading reports for the commit 611a124 to get more accurate results

Files Patch % Lines
src/Mechanical/PlanarMechanics/sensors.jl 0.00% 263 Missing :warning:
src/Mechanical/PlanarMechanics/joints.jl 0.00% 45 Missing :warning:
src/Mechanical/PlanarMechanics/components.jl 0.00% 31 Missing :warning:
Additional details and impacted files
@@             Coverage Diff             @@
##             main     #220       +/-   ##
===========================================
- Coverage   54.12%   12.98%   -41.15%     
===========================================
  Files          48       52        +4     
  Lines        1648     1987      +339     
===========================================
- Hits          892      258      -634     
- Misses        756     1729      +973     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

codecov[bot] avatar Sep 23 '23 06:09 codecov[bot]

I think this idea is great! Down the line we may move this out to a separate SciML/PlanarMechanics to reduce the loading time of the standard library, but we can do any reformatting like that down the line so I wouldn't worry about that right now. For now, just getting more modules is good and we'd be happy to help maintain them.

ChrisRackauckas avatar Sep 24 '23 14:09 ChrisRackauckas

Is there a way I can do this in MTK, i.e., asserting the number of connections?

assert(cardinality(frame_a) > 0, "Connector frame_a must be connected at least once");

AbdulrhmnGhanem avatar Sep 28 '23 07:09 AbdulrhmnGhanem

Is there a way I can do this in MTK, i.e., asserting the number of connections?

Not at the moment, but it would be a nice feature to have

baggepinnen avatar Nov 08 '23 05:11 baggepinnen

Not at the moment, but it would be a nice feature to have

Could you hint where I can start to add this in MTK?

AbdulrhmnGhanem avatar Nov 17 '23 23:11 AbdulrhmnGhanem

Could you hint where I can start to add this in MTK?

I'm not sure TBH, but this feature is a debugging feature, the PR here should be good to go without it.

Do you feel that the PR has reached a state where it's ready for review?

baggepinnen avatar Nov 20 '23 03:11 baggepinnen

Do you feel that the PR has reached a state where it's ready for review?

Sorrowfully no, but I will resume working on it and get it ready as soon as possible.

AbdulrhmnGhanem avatar Nov 20 '23 14:11 AbdulrhmnGhanem

@baggepinnen, @YingboMa, could you take a look at the pendulum model below? https://github.com/SciML/ModelingToolkitStandardLibrary.jl/blob/7051dbea205126115cd56e21850c30eb0ee17c83/test/Mechanical/planar_mechanics.jl#L35-L61

I have been stuck with it for quite long time so far.

Why would such a system have 7 states after simplification? it's unlikely to be a problem with structural_simplify but I ran out of debugging tricks.

AbdulrhmnGhanem avatar Dec 20 '23 14:12 AbdulrhmnGhanem

Why would such a system have 7 states after simplification?

What are these 7 variables?

These equations in Revolute look suspicious

phi ~ ifelse(constant_phi === nothing, phi, constant_phi)

You are adding the equation phi ~ phi here, I'm not sure if that trips MTK up?

The constant_phi argument is not documented, so I'm not sure what it's supposed to accomplish. If you keep the angle constant, it's not really a joint any more?

baggepinnen avatar Mar 07 '24 09:03 baggepinnen

What are these 7 variables?

julia> states(sys)
7-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 rod₊frame_b₊y(t)
 rod₊frame_b₊yˍt(t)
 rod₊frame_a₊phi(t)
 rod₊frame_b₊fx(t)
 rod₊frame_a₊phiˍt(t)
 rod₊frame_a₊phiˍtt(t)
 body₊ryˍtt(t)

julia> equations(sys)
7-element Vector{Equation}:
 Differential(t)(rod₊frame_b₊y(t)) ~ rod₊frame_b₊yˍt(t)
 Differential(t)(rod₊frame_b₊yˍt(t)) ~ body₊ryˍtt(t)
 0 ~ rod₊frame_b₊y(t) - rod₊frame_a₊y(t) - rod₊rx*sin(rod₊frame_a₊phi(t)) - rod₊ry*cos(rod₊frame_a₊phi(t))
 0 ~ rod₊frame_b₊j(t) + body₊frame₊j(t)
 0 ~ rod₊frame_b₊xˍtt(t) - rod₊frame_a₊xˍtt(t) + rod₊rx*sin(rod₊frame_a₊phi(t))*rod₊frame_a₊phiˍtt(t) + rod₊ry*rod₊frame_a₊phiˍtt(t)*cos(rod₊frame_a₊phi(t)) + rod₊rx*cos(rod₊frame_a₊phi(t))*(rod₊frame_a₊phiˍt(t)^2) - rod₊ry*sin(rod₊frame_a₊phi(t))*(rod₊frame_a₊phiˍt(t)^2)
 0 ~ rod₊frame_b₊yˍt(t) - rod₊frame_a₊yˍt(t) - rod₊rx*cos(rod₊frame_a₊phi(t))*rod₊frame_a₊phiˍt(t) + rod₊ry*sin(rod₊frame_a₊phi(t))*rod₊frame_a₊phiˍt(t)
 0 ~ -rod₊frame_a₊yˍtt(t) + body₊ryˍtt(t) - rod₊rx*rod₊frame_a₊phiˍtt(t)*cos(rod₊frame_a₊phi(t)) + rod₊ry*sin(rod₊frame_a₊phi(t))*rod₊frame_a₊phiˍtt(t) + rod₊rx*sin(rod₊frame_a₊phi(t))*(rod₊frame_a₊phiˍt(t)^2) + rod₊ry*cos(rod₊frame_a₊phi(t))*(rod₊frame_a₊phiˍt(t)^2)

You are adding the equation phi ~ phi here, I'm not sure if that trips MTK up?

Removing it didn't make a difference.

The constant_phi argument is not documented, so I'm not sure what it's supposed to accomplish. If you keep the angle constant, it's not really a joint any more?

You're right, I needed constant_ω to make a revolute move with a constant angular velocity and added constant_phi, and constat_tau just in case. Is there a better way to do this?

https://github.com/SciML/ModelingToolkitStandardLibrary.jl/blob/abc9d7f7334ee9c151f03f7970ff1ca04a2c600d/test/Mechanical/planar_mechanics.jl#L79

AbdulrhmnGhanem avatar Mar 07 '24 09:03 AbdulrhmnGhanem

It looks like MTK has chosen rod variables as state

 rod₊frame_a₊phi(t)
 rod₊frame_a₊phiˍt(t)

these are equivalent to the revolute joint angle and angular velocity and should be sufficient as state variables.

You're right, I needed constant_ω to make a revolute move with a constant angular velocity and added constant_phi, and constat_tau just in case. Is there a better way to do this?

Something like this component?

https://github.com/modelica/ModelicaStandardLibrary/blob/master/Modelica/Mechanics/Rotational/Sources/ConstantSpeed.mo

i.e., a constant-speed source

baggepinnen avatar Mar 07 '24 11:03 baggepinnen

The problem here might be the same problem that we have with the multibody library. Multibody requires JuliaSimCompiler to simpify properly. I tried the pendulum example with JuliaSimCompiler but ran into a bug, so I cannot verify that this is indeed the case

julia> sys = structural_simplify(JuliaSimCompiler.IRSystem(model))
ERROR: BoundsError: attempt to access 61-element Vector{Vector{Int64}} at index [62]
Stacktrace:
 [1] getindex
   @ ./essentials.jl:13 [inlined]
 [2] 𝑑neighbors

baggepinnen avatar Mar 07 '24 13:03 baggepinnen

i.e., a constant-speed source

Yes, this is better, I'll add it.

AbdulrhmnGhanem avatar Mar 07 '24 14:03 AbdulrhmnGhanem

The problem here might be the same problem that we have with the multibody library.

I merged the main branch into this branch, just in case this bug was fixed in the latest version.

AbdulrhmnGhanem avatar Mar 08 '24 13:03 AbdulrhmnGhanem