Lumped beam model

Hello,
I’m trying to implement a 2D lumped beam model with modelingtoolkit. Using the small angle approximation, I managed to have something working with 5 elements. But increasing the number of elements to 6 make ODEProblem last forever …any clues ?

Plus when I try to use exact version with cos and sine (you can try by uncommenting the other rigidlink definition) I received the following error msg:

ERROR: MethodError: no method matching var_derivative!(::TearingState{NonlinearSystem}, ::Int64)

Any ideas on the two topics greatly appreciated !

using ModelingToolkit
using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkit: t_nounits as t, D_nounits as D
import DifferentialEquations
using Plots
#-----------------------------------custom 2D mech library----------------------------------
@connector Frame begin
    x(t), [description = "x position of frame"]
    y(t), [description = "y position of frame"]
    phi(t), [description = "Rotation angle of frame"]
    tau(t), [connect=Flow, description = "Cut torque in frame"]
    fx(t), [connect=Flow, description = "Cut force x in frame"]
    fy(t), [connect=Flow, description = "Cut force y in frame"]
end

@mtkmodel Fixed begin
    @components begin
        frame = Frame()
    end
    @parameters begin
        x0 = 0.0, [description = "Fixed x position"]
        y0 = 0.0, [description = "Fixed y position"]
        phi0 = 0.0, [description = "Fixed y position"]
    end
    @equations begin
        frame.x ~ x0
        frame.y ~ y0
        frame.phi ~ phi0
    end
end

@mtkmodel PartialCompliantRotational begin
    @components begin
        frame_a = Frame()
        frame_b = Frame()
    end
    @variables begin
        phi_rel(t), [description = "Relative rotation angle between frames", guess = 0.0]
        w_rel(t), [description = "Relative angular velocity between flanges", guess = 0.0]
        tau(t), [description = "Torque between frames", guess = 0.0]
    end
    @equations begin
        phi_rel ~ frame_b.phi - frame_a.phi
        D(phi_rel) ~ w_rel
        frame_b.tau ~ tau
        frame_a.tau ~ -tau
        frame_a.x ~ frame_b.x
        frame_a.y ~ frame_b.y
        frame_a.fx ~ -frame_b.fx
        frame_a.fy ~ -frame_b.fy
    end
end

@mtkmodel SpringDamperRotational begin
    @extend phi_rel, w_rel, tau = partial_comp = PartialCompliantRotational()
    @variables begin
        tau_c(t), [description = "Spring torque"]
        tau_d(t), [description = "Damper torque"]
    end
    @parameters begin
        d, [description = "Damping constant"]
        c, [description = "Spring constant"]
        phi_rel0 = 0.0, [description = "Unstretched spring angle"]
    end
    @equations begin
        tau_c ~ c * (phi_rel - phi_rel0)
        tau_d ~ d * w_rel
        tau ~ tau_c + tau_d
    end
end

@mtkmodel MassInertia begin
    @components begin
        frame = Frame()
    end
    @parameters begin
        m = 1.0, [description = "Mass kg"]
        i = 10.0, [description = "Inertia kg⋅m²"]
    end
    @variables begin
        x(t), [description = "Absolute x position", guess = 0.0]
        vx(t), [description = "Absolute x velocity", guess = 0.0]
        ax(t), [description = "Absolute x acceleration", guess = 0.0]
        y(t), [description = "Absolute y position", guess = 0.0]
        vy(t), [description = "Absolute y velocity", guess = 0.0]
        ay(t), [description = "Absolute y acceleration", guess = 0.0]
        phi(t), [description = "Absolute angular position", guess = 0.0]
        w(t), [description = "Absolute angular velocity", guess = 0.0]
        a(t), [description = "Absolute angular acceleration", guess = 0.0]
    end
    @equations begin
        frame.x ~ x
        D(x) ~ vx
        D(vx) ~ ax
        frame.y ~ y
        D(y) ~ vy
        D(vy) ~ ay
        frame.phi ~ phi
        D(phi) ~ w
        D(w) ~ a
        m*ax ~ frame.fx 
        m*ay ~ frame.fy 
        i*a ~ frame.tau 
    end
end

# @mtkmodel RigidLink begin
#     @parameters begin
#         l = l
#         phi_0 = 0, [description = "link angle"]
#     end
#     @components begin
#         frame_a = Frame()
#         frame_b = Frame()
#     end
#     @equations begin
#         # geometric
#         frame_b.x ~ frame_a.x + l*cos(phi_0+frame_a.phi)
#         frame_b.y ~ frame_a.y + l*sin(phi_0+frame_a.phi)
#         frame_b.phi ~ frame_a.phi+phi_0
#         # torque
#         frame_a.tau ~ -(-frame_b.fx*l*sin(phi_0+frame_a.phi) + 
#             frame_b.fy*l*cos(phi_0+frame_a.phi)) - frame_b.tau
#         # forces
#         frame_a.fx ~ -frame_b.fx
#         frame_a.fy ~ -frame_b.fy
#     end
# end

@mtkmodel RigidLink begin
    @parameters begin
        l = l
    end
    @components begin
        frame_a = Frame()
        frame_b = Frame()
    end
    @equations begin
        # geometric
        frame_b.x ~ frame_a.x + l*(1)
        frame_b.y ~ frame_a.y + l*(frame_a.phi)
        frame_b.phi ~ frame_a.phi
        # torque
        frame_a.tau ~ -(-frame_b.fx*l*(frame_a.phi) + 
            frame_b.fy*l*(1)) - frame_b.tau
        # forces
        frame_a.fx ~ -frame_b.fx
        frame_a.fy ~ -frame_b.fy
    end
end

@mtkmodel VariableForce begin
    @components begin
        fx = RealInput()
        fy = RealInput()
        frame = Frame()
    end
    @equations begin
        frame.fx ~ -fx.u
        frame.fy ~ -fy.u
        frame.tau ~ 0.0
    end
end

@mtkmodel BeamElement begin
    @parameters begin
        L = 50e-3, [description = "length of element"] # m
        B = 0.8e-3, [description = "width of beam"] # m
        H = 20.9e-3, [description = "height of beam"] # m
        E = 200.0e9 , [description = "Modulus"] 
        ρ = 8000.0, [description = "Density"]
        d = 0.01, [description = "damping"]
    end
    @components begin
        spring = SpringDamperRotational(name=:spring,phi_rel0=0.0,c=E*(H*B^3/12)/L,d=d)
        beam_1 = RigidLink(;name=:beam_1,l=L/2)
        beam_2 = RigidLink(;name=:beam_2,l=L/2)
        mass = MassInertia(;name=:mass,m=L*B*H*ρ,i=(L*B*H*ρ)*(L^2+B^2)/12)
    end
    @equations begin
        connect(spring.frame_b,beam_1.frame_a)
        connect(beam_1.frame_b,beam_2.frame_a,mass.frame) #mass.frame
    end
end

#-------------------------------System definition------------------------------------------#
N=5
@mtkmodel system begin
    @components begin
        fixed = Fixed(;name=:fixed)
        elements = [BeamElement(;name=Symbol(string("el_",i))) for i in 1:N]
        force = VariableForce(;name=:force)
        step = Step(;name=:step, height=1, offset=0, start_time=0.1, duration=Inf, smooth=true)
    end
    @equations begin
        # for i=1:1 # use below syntax
        #     connect(elements[i].beam_2.frame_b,elements[i+1].spring.frame_a)    
        # end
        [connect(elements[i].beam_2.frame_b,elements[i+1].spring.frame_a) for i in 1:(N-1)]
        connect(fixed.frame,elements[1].spring.frame_a)
        connect(elements[N].beam_2.frame_b,force.frame)
        connect(step.output,force.fy)
        force.fx.u ~ 0.0
    end
end

@mtkbuild sys = system()
@time prob = ODEProblem(sys,[[eval(Meta.parse("sys.el_$i.spring.phi_rel => 0.0")) for i in 1:N]
    [eval(Meta.parse("sys.el_$i.spring.w_rel => 0.0")) for i in 1:N];]
    , (0, 1.0)) #[eval(Meta.parse("sys.el_$i.spring.w_rel => 0.0")) for i in 1:3]]

sol = solve(prob)
plot(sol,idxs=[eval(Meta.parse("sys.el_$i.beam_2.frame_b.y")) for i in 1:N])

My guess is that you get a poorly balanced statespace realization that causes the solvers to struggle. You could verify whether or not this is the case by running

@named model = system()
mats, ssys = linearize(model, [], [])
mats.A

and inspect this A matrix. It could also be that MTK fails to simplify the system properly and you get a large-dimensional DAE instead of the ODE you’d expect, this does indeed appear to be the case


julia> equations(sys)
21-element Vector{Equation}:
 Differential(t)(el_5₊mass₊frame₊phi(t)) ~ el_5₊mass₊frame₊phiˍt(t)
 Differential(t)(el_4₊mass₊frame₊phi(t)) ~ el_4₊mass₊frame₊phiˍt(t)
 Differential(t)(el_3₊mass₊frame₊phi(t)) ~ el_3₊mass₊frame₊phiˍt(t)
 Differential(t)(el_3₊spring₊frame_a₊phi(t)) ~ el_3₊spring₊frame_a₊phiˍt(t)
 Differential(t)(el_5₊mass₊frame₊y(t)) ~ el_5₊mass₊frame₊yˍt(t)
 Differential(t)(el_3₊spring₊frame_a₊phiˍt(t)) ~ el_2₊mass₊phiˍtt(t)
 Differential(t)(el_3₊mass₊frame₊phiˍt(t)) ~ el_3₊mass₊phiˍtt(t)
 Differential(t)(el_4₊mass₊frame₊phiˍt(t)) ~ el_4₊beam_2₊frame_a₊phiˍtt(t)
 Differential(t)(el_5₊mass₊frame₊phiˍt(t)) ~ el_5₊beam_2₊frame_a₊phiˍtt(t)
 Differential(t)(el_5₊mass₊frame₊yˍt(t)) ~ el_5₊mass₊yˍtt(t)
 0 ~ -el_2₊spring₊frame_a₊phi(t) + el_1₊beam_2₊frame_a₊phi(t)
 0 ~ -el_5₊mass₊frame₊fy(t) - el_5₊beam_1₊frame_b₊fy(t) + step₊output₊u(t)
 0 ~ el_5₊mass₊frame₊fy(t) - el_5₊mass₊m*el_5₊mass₊ay(t)
 0 ~ -el_1₊mass₊frame₊phiˍt(t) + el_1₊mass₊phiˍt(t)
 0 ~ el_1₊mass₊frame₊phiˍtt(t) - el_1₊beam_1₊frame_a₊phiˍtt(t)
 0 ~ -el_2₊mass₊phiˍtt(t) + el_2₊beam_2₊frame_a₊phiˍtt(t)
 0 ~ -el_3₊beam_1₊frame_b₊yˍtt(t) + el_3₊beam_2₊frame_a₊yˍtt(t)
 0 ~ -el_3₊mass₊phiˍtt(t) + el_3₊beam_2₊frame_a₊phiˍtt(t)
 0 ~ el_4₊mass₊phiˍtt(t) - el_4₊beam_2₊frame_a₊phiˍtt(t)
 0 ~ el_5₊beam_2₊frame_a₊yˍtt(t) - el_5₊mass₊yˍtt(t)
 0 ~ -el_5₊beam_2₊frame_a₊phiˍtt(t) + el_5₊mass₊phiˍtt(t)

Have you tried specifying a different solver manually?

Hello,
Actually the solving of the system is ok, it is the “construction/initialisation” phase i.e. the ODEProblem instruction that blocks for more than 5 elements:

prob = ODEProblem(sys,[[eval(Meta.parse("sys.el_$i.spring.phi_rel => 0.0")) for i in 1:N]
    [eval(Meta.parse("sys.el_$i.spring.w_rel => 0.0")) for i in 1:N];]
    , (0, 1.0))

I tried what you suggest but it seems not able to linearize the problem, the last error msg repeat indefinitely…:

mats, ssys = linearize(model, [],[])
┌ Warning: An empty operating point was passed to `linearization_function`. An operating point containing the variables that will be changed in `linearize` should be provided. Disable this warning by passing `warn_empty_op = false`.
└ @ ModelingToolkit C:\Users\yjr\.julia\packages\ModelingToolkit\aau6A\src\linearization.jl:51
┌ Warning: The initialization system is structurally singular. Guess values may significantly affect the initial values of the ODE. The problematic variables are Any[el_5₊beam_2₊frame_a₊yˍt(t), force₊frame₊y(t), el_2₊beam_1₊frame_b₊y(t), el_2₊spring₊frame_a₊tau(t), el_1₊beam_1₊frame_b₊tau(t), el_5₊mass₊a(t), el_4₊mass₊a(t), el_2₊beam_1₊frame_b₊tau(t), el_2₊beam_2₊frame_a₊fy(t), el_3₊beam_1₊frame_a₊phiˍtt(t)].
│ 
│ Note that the identification of problematic variables is a best-effort heuristic.
└ @ ModelingToolkit C:\Users\yjr\.julia\packages\ModelingToolkit\aau6A\src\systems\diffeqs\abstractodesystem.jl:1493
┌ Warning: Initialization system is underdetermined. 9 equations for 19 unknowns. Initialization will default to using least squares. `SCCNonlinearProblem` can only be used for initialization of fully determined systems and hence will not be used here. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true
│ 
│ Note that the identification of problematic variables is a best-effort heuristic.
└ @ ModelingToolkit C:\Users\yjr\.julia\packages\ModelingToolkit\aau6A\src\systems\diffeqs\abstractodesystem.jl:1493
│ 
│ Note that the identification of problematic variables is a best-effort heuristic.
│ 
│ 
│ Note that the identification of problematic variables is a best-effort heuristic.
└ @ ModelingToolkit C:\Users\yjr\.julia\packages\ModelingToolkit\aau6A\src\systems\diffeqs\abstractodesystem.jl:1493
┌ Warning: Initialization system is underdetermined. 9 equations for 19 unknowns. Initialization will default to using least squares. `SCCNonlinearProblem` can only be used for initialization of fully determined systems and hence will not be used here. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true
┌ Warning: Initialization system is underdetermined. 9 equations for 19 unknowns. Initialization will default to using least squares. `SCCNonlinearProblem` can only be used for initialization of fully determined systems and hence will not be used here. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true
e will not be used here. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true
└ @ ModelingToolkit C:\Users\yjr\.julia\packages\ModelingToolkit\aau6A\src\systems\diffeqs\abstractodesystem.jl:1523
┌ Warning: Did not converge after `maxiters = 100` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics C:\Users\yjr\.julia\packages\Symbolics\tSVcV\src\variable.jl:587
┌ Warning: Did not converge after `maxiters = 100` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.

I have a similar issue when modeling a 1D flow channel. I have 100 elements and connecting the inlet to the outlet exactly as you have done. Then ODEproblem takes 3 minutes, which is not OK. I found the ‘solution’ is to make use of arrays. Instead of having 100 tau_c(t) states, you have one (tau_c(t))[1:100] array state. The only issue is that you now need to model equation-based and not component-based, which sort of defeats the purpose of using MTK. For me, using arrays, ODEProblem went from 3 minutes to 10 seconds.

I am happy to share my code if that might help?

For the sin issue, I don’t know exactly what MTK understands to be symbolic and what not, but it could help to @register_symbolic sin so that MTK does not try to undertsand what sin means, but sees it as a black-box sort of.

that’s already done by default

1 Like

The problem may have something to do with the modelling of the mechanical system.
I tried to run the example with MODIA and encountered a similar problem. N=4 ran smoothly, N=5 resulted in an error message (below).

Error message from getSortedAndSolvedAST for model Lumped:
Cannot transform to ODE because constraint equations cannot be explicitly solved by tearing
Involved variables:
    element1.mass.phi
    element2.mass.phi
    element5.mass.phi
...
...

With a modification of MassInertia I was able to run your example.

@mtkmodel MassInertia begin
    @components begin
        frame = Frame()
    end
    @parameters begin
        m = 1.0, [description = "Mass kg"]
        i = 10.0, [description = "Inertia kg⋅m²"]
    end
    @variables begin
        x(t), [description = "Absolute x position", guess = 0.0]
        vx(t), [description = "Absolute x velocity", guess = 0.0]
        ax(t), [description = "Absolute x acceleration", guess = 0.0]
        y(t), [description = "Absolute y position", guess = 0.0]
        vy(t), [description = "Absolute y velocity", guess = 0.0]
        ay(t), [description = "Absolute y acceleration", guess = 0.0]
        phi(t), [description = "Absolute angular position", guess = 0.0]
        w(t), [description = "Absolute angular velocity", guess = 0.0]
#         a(t), [description = "Absolute angular acceleration", guess = 0.0]
    end
    @equations begin
        frame.x ~ x
        D(x) ~ vx
        D(vx) ~ ax
        frame.y ~ y
        D(y) ~ vy
        D(vy) ~ ay
        frame.phi ~ phi
        D(phi) ~ w
#         D(w) ~ a
        m*ax ~ frame.fx
        m*ay ~ frame.fy
#         i*a ~ frame.tau
        0.0 ~ frame.tau
    end
end

This is not the solution, but it can be helpful in finding the main problem.