Convert string to equation (fast)

General Usage
,
1

Hi, I’m pretty new to the Julia language and i have a problem with my program.

Just a short background, I am a engineering student and my program shall be a Multibody simulation tool. I have already solved the mechanics and necessary equations with the program Mathematica.
The goal is to put these equations in a loop (in Julia) in which i calculate some variables and calculate the equations for the next step to get the trajectiories of my Bodies.
These equations are in a .txt file (they are pretty long, ~12 pages for 3 equations in pdf) that i read at the beginning of the program as a string.
I also want to be flexible if some boundary conditions change, so i dont want to hardcode the equations.

I tried to parse the equations with eval(Meta.parse(“equation”)) outside of my loop but then in the loop i have no possibility to recalculate this equation with the uptdated variables. When i try to evaluate everything inside the loop it takes way too long for my simulation (for 1600 steps it takes ~4min → my program shall be real time).
Is there another way to “save” my equations beforehand so that i can evaluate everything much quicker?

Im not sure if i shall put this question in the “Performance” category because I think there has to be a simple solution which i havent found yet.

Thanks in advance,
Daniel

2

Hi Daniel,

Welcome to the Julia community!

Could you condense your current approach into a minimal working example (with simplified equations and update rules)? (See also this PSA.)
Currently it is not clear to me why you could not just load your equations once as functions of the updating variables and then simply continually evaluate these functions in your loop.

3

Thanks, of course, i should have added the structure in the first post…
This is the part where i read the equations and part of my simulation loop. This is not the whole loop, it “works” but the f_dDG2/4/6 do not get udated in every loop, just in the first one when i parse it. I hope it is understandable…
Thanks!

## Equations 
eq_motion1 = readlines(path_eq_motion)
eq_motion = split(eq_motion1[1], ",")
# Equations
f_dDG2 = eval(Meta.parse(eq_motion[1]))
f_dDG4 = eval(Meta.parse(eq_motion[2]))
f_dDG6 = eval(Meta.parse(eq_motion[3]))

# f_dDG2:  "DG2 = dDG2 = -((-((-((2*l3*l5*m3*cos(DG3)*cos(DG5)..."
# The equations are very long and consist mostly of 
# constants (calculated in step t-1) and trigonimic functions


## SIMULATION LOOP
for step in 1:n
    println(step)
    # Calculate the current time
    t = t0 + step * Δt
    # Access previous values from results df or initial conditions
    if step > 1
        x1_p = results.DG1[end]         # x1    
        dx1_p = results.DG2[end]        # x1'
        φ_p = results.DG3[end]          # φ
        dφ_p = results.DG4[end]         # φ'
        ψ_p = results.DG5[end]          # ψ
        dψ_p = results.DG6[end]         # ψ'
    else
        # initial conditions
        x1_p = 0                        # x1 at t=0
        dx1_p = -56.0/3.6          # x1' at t=0
        φ_p = deg2rad(20.0)             # φ at t=0
        dφ_p = 0.0                      # φ' at t=0
        ψ_p = 0.0                       # ψ at t=0
        dψ_p = 0.0                      # ψ' at t=0
    end

    # check the boundaries of Input acceleration
        if t < last(acc.Time)
            ax4 = lin_interpolate(acc, :Time, :Acc, t)
        else
            println("Simulation finished -> no acc data for that timestep")
            break
        end
    

    # Euler method: Calculate new states
    # Calculate accelerations
    dDG2 = f_DG2  # function for x1''
    dDG4 = f_DG4  # function for φ''
    dDG6 = f_DG6  # function for ψ''
    # Update velocities and positions (Euler method)
    global DG2 = DG2 + Δt * dDG2;  # Velocity of x1
    global DG1 = DG1 + Δt * DG2;   # Position of x1 
    global DG4 = DG4 + Δt * dDG4;  # Velocity of φ
    global DG3 = DG3 + Δt * DG4;   # Position of φ
    global DG6 = DG6 + Δt * dDG6;  # Velocity of ψ
    global DG5 = DG5 + Δt * DG6;   # Position of ψ
    # Store current values in the DataFrame
    push!(results, (step, t, ax4, DG1, DG2, DG3, DG4, DG5, DG6));
end
4

Why must the equations be written in a text-file instead of as program code in a jl-file?

It’s still hard-coded, only in txt format.

Are you updating the files during program execution? How, and what sort of changes? Can’t this update be part of the program?

2 Likes
5

This is part of my bachelor thesis and i could not figure out a way to restructure and solve these equations in julia. So i switched to mathematica (where i have some experience) to solve the equations. In the long term (this is part of a bigger project) it shall be standalone in julia but for now i do not have enough time.
I am currently working out all the small details and validating the equations. So i need a “quick” way to check if my trajectories work out.

I dont think that there is a way to export the outputs from mathematica directly to julia.

With hardcode i meant the equations written in the julia file, with the txt file i can make some changes in the mathematica file and update the output txt.

I do not have a lot of experience with coding so i know there are better ways to code/execute so i am oped for sugestions.

During the execution of the julia program, i do not need to change the equations, just in general.

Thanks!

6

But if you can convert the txt file to Julia code simply with split and eval(Meta.parse()), then it must almost already be valid Julia, right? It should be equally straightforward to turn it into a jl-file once, as parsing it repeatedly?

Could you share a snippet of this txt file, or a dummy version?

7

Another question is whether it’s even just the parsing that’s causing the slowness in the first place. The code snippet you shared is apparently using lots of globals, and not wrapping your main loop in a function. This will be very slow in any case.

Take a look at the performance tips section of the manual: Performance Tips · The Julia Language

2 Likes
8

Thanks for the code snippet.

I’m afraid it’s not a proper MWE though, so some unclarity remains.

Should we interpret this as the String eq_motion[1] being "DG2 = dDG2 = <...>"? If so, where and how did you define all the variables in the right-hand side (l3, l5, m3, DG3, etc.)? Also, note that f_dDG2 (assuming it is indeed eval(Meta.parse(<that string>))) evaluates to the (evaluated) right-hand side, but I would consider this a side-effect. Indeed, you already assign a value to DG2 and dDG2, so that the dDG2 = f_DG2 (where I’m going to pretend f_DG2 == f_dDG2, as you did not include a declaration) probably is not needed in step==1.

E.g.

julia> f = eval(Meta.parse("x = y = 2"))
2

julia> @show x y f  # Note that x and y are defined. No need for x = f
x = 2
y = 2
f = 2
2


The function idea I suggested was to create a function f_DG2 once, via something like
f_DG2 = eval(Meta.parse("(l3, l5, DG3, DG5, <...>) -> -((-((-((2*l3*l5*<...>")
once, and then in your loop update DG2 as
DG2 = f_DG2(l3, l5, DG3, DG5, <...>)
using the previously computed values of the variables. E.g.

julia>  f_x = eval(Meta.parse("(a, b) -> (a + 1)^2 - b"))
#1 (generic function with 1 method)

julia> for (a, b) in zip(1:3, 10:10:30) 
           x = f_x(a, b); @show a b x; println()
       end
a = 1
b = 10
x = -6

a = 2
b = 20
x = -11

a = 3
b = 30
x = -14
9

Hi,

yeah i didnt really understand the “MWE”…
Down below i added the now (hopefully) working MWE.

Yes i tried to do so…
In the example below the equations are similar to my real equations, just a lot shorter.

This code below works like my original one, and the dDG2/4/6 do not get updated:

using DataFrames
using Interpolations
using Plots

function lin_interpolate(df::DataFrame, col1::Symbol, col2::Symbol, col1_value::Float64)
    # extract the columns in which you want to interpolate
    x = collect(df[!, col1])
    y = collect(df[!, col2])

    # Check if the col1_value is within the range of the dataframe
    if col1_value < minimum(x) || col1_value > maximum(x)
        return "error"
    end
    # linear interpolate (Module)
    itp = linear_interpolation(x, y)
    # calculate the interpolation for var1_value
    col2_value = itp(col1_value)
    # return the interpolated value
    return col2_value
end


## BODIES - GEOMETRICAL VALUES
    # Initialize the geometrical values
    # First Body
    l1 = 20.0
    l2 = 20.0
    m1 = 20.0
    φ0 = deg2rad(20.0)
    ψ0 = deg2rad(0.0)

    # Initial length of the springs
    x10 = 0.0
    l0Fbt = sqrt((l1+x10*sin(φ0))^2+(l2+x10*cos(φ0)))
    J10 = φ0
    J20 = deg2rad(90.0) + φ0 - ψ0
# 
## FORCES - BELT/AIRBAG
    # Initialise lookup tables for the forces
    Fbt_v = DataFrame(Distance = Float64[], Force = Float64[])
    push!(Fbt_v, (-100.0, 0.0))
    push!(Fbt_v, (0.0, 0.0))
    push!(Fbt_v, (100.0, 0.0))
    push!(Fbt_v, (200.0, 3500.0))
    push!(Fbt_v, (1000.0, 3500.0))
#
## TORQUES - JOINTS
    # Initialise lookup tables for the torques
    Mg1_v = DataFrame(Angle = Float64[], Torque = Float64[])
    push!(Mg1_v, (deg2rad(-180.0), 20.0))
    push!(Mg1_v, (deg2rad(180.0), 20.0))
#
## SIMULATION VARIABLES
    # save the input Acceleration in a DataFrame
    acc = DataFrame(Time = Float64[], Acc = Float64[])
    push!(acc, (0.0, 0.0))
    push!(acc, (0.05, 200.0))

    Δt = 2.5e-5                                
    T = round.(last(acc.Time), digits=5)        # Length of Simulation  
    n = T / Δt                                  # Number of Steps the Simulation takes -> calculated from input acc data
    t0 = round.(acc[1,1], digits=5)             # Starting time

    results = DataFrame(Steps=Int64[], Time = Float64[], Acc= Float64[], DG1 = Float64[], DG2 = Float64[], DG3 = Float64[], DG4 = Float64[], DG5 = Float64[], DG6 = Float64[], dDG2 = Float64[], dDG4 = Float64[], dDG6 = Float64[], dFbt = Float64[], Fbt = Float64[], dJ1 = Float64[], Mg1 = Float64[]);
#
## EQUATIONS 
    ## Initialize the equations
    eq_motion1 = "dDG2=2*l1*m1*cos(DG3)*cos(DG5)*DG1+DG4*DG6+Mg1, dDG4=l1*m1*cos(DG3)*tan(DG5)*DG4-l1*m1*cos(DG3)*cos(DG5)*DG6+Mg1, dDG6=2*l1*m1*cos(DG3)*sin(DG5)*DG4+DG2/DG6+Mg1"
    eq_motion = split(eq_motion1, ",")

    # Equations string
    f_dDG2 = eq_motion[1]
    f_dDG4 = eq_motion[2]
    f_dDG6 = eq_motion[3]

    parsed_dDG2 = Meta.parse(f_dDG2)
    parsed_dDG4 = Meta.parse(f_dDG4)
    parsed_dDG6 = Meta.parse(f_dDG6)

    # Equations
    c_dDG2 = eval(parsed_dDG2)
    c_dDG4 = eval(parsed_dDG4)
    c_dDG6 = eval(parsed_dDG6)
#
## SIMULATION LOOP
    for step in 1:n
        println(step)
        # Calculate the current time
        t = t0 + step * Δt
        # Access previous values from results DataFrame if we are not at the first step
        if step > 1
            DG1 = results.DG1[end]          
            DG2 = results.DG2[end] 
            DG3 = results.DG3[end] 
            DG4 = results.DG4[end]         
            DG5 = results.DG5[end]           
            DG6 = results.DG6[end]  
            x1_p = results.DG1[end]         # x1    
            dx1_p = results.DG2[end]        # x1'
            φ_p = results.DG3[end]          # φ
            dφ_p = results.DG4[end]         # φ'
            ψ_p = results.DG5[end]          # ψ
            dψ_p = results.DG6[end]         # ψ'
        else
            # Use initial conditions for the first step
            DG1 = 0.0;          
            DG2 = 56.0*1000/3.6;
            dDG2 = 0.0
            DG3 = deg2rad(20.0);
            DG4 = 0.0;    
            dDG4 = 0.0;      
            DG5 = 0.0;          
            DG6 = 0.0; 
            dDG6 = 0.0;         
            x1_p = DG1;                                 # x1 at t=0
            dx1_p = DG2;                     #   x1' at t=0
            φ_p = DG3;                       # φ at t=0
            dφ_p = DG4;                                 # φ' at t=0
            ψ_p = DG5;                                   # ψ at t=0
            dψ_p = DG6;                                 # ψ' at t=0
        end
        # Calculate the lenghts of the springs/angles 
        lFbt = sqrt((l1+x1_p*sin(φ_p))^2+(l2+x1_p*cos(φ_p)))
        J1 = φ_p
        # Calculate the length difference of the springs to t0
        dlFbt = lFbt - l0Fbt
        dJ1 = J1 - J10
        # Calculate the t-1 forces/Torques
        Fbt = lin_interpolate(Fbt_v, :Distance, :Force, dlFbt)
        Mg1 = lin_interpolate(Mg1_v, :Angle, :Torque, dJ1)
    
        # check the boundaries of Input acc
            if t < last(acc.Time)
                ax4 = lin_interpolate(acc, :Time, :Acc, t)
                #println(acc_a)
            else
                println("Simulation finished -> no acc data for that timestep")
                break
            end
        #
        # calculate the Forces
        # Euler method: Calculate new states
        # Calculate accelerations
        dDG2 = c_dDG2
        dDG4 = c_dDG4
        dDG6 = c_dDG6
        # Update velocities and positions (Euler method)
        DG2 = DG2 + Δt * dDG2;  # Velocity of x1
        DG1 = DG1 + Δt * DG2;   # Position of x1 
        DG4 = DG4 + Δt * dDG4;  # Velocity of φ
        DG3 = DG3 + Δt * DG4;   # Position of φ
        DG6 = DG6 + Δt * dDG6;  # Velocity of ψ
        DG5 = DG5 + Δt * DG6;   # Position of ψ
        # Store current values in the DataFrame
        push!(results, (step, t, ax4, DG1, DG2, DG3, DG4, DG5, DG6, dDG2, dDG4, dDG6, dlFbt, Fbt, dJ1, Mg1));
    end

    rename!(results,
    :DG1 => :x1,
    :DG2 => :vx1,
    :DG3 => :Phi,
    :DG4 => :vPhi,
    :DG5 => :Psi,
    :DG6 => :vPsi
    )

    ## PLOTS
    plot(results.Time, results.Acc, label="acc", xlabel="Time (s)", ylabel="", legend=:topright)
    plot!(results.Time, results.x1, label="x1", linestyle=:dash)
    plot!(results.Time, results.vx1, label="vx1", linestyle=:dash)
    plot!(results.Time, results.Phi, label="Phi", linestyle=:dash)
    plot!(results.Time, results.vPhi, label="vPhi", linestyle=:dash)
    plot!(results.Time, results.Psi, label="Psi", linestyle=:dash)
    plot!(results.Time, results.vPsi, label="vPsi", linestyle=:dash)
10

I think it is the parsing/evaluating…

I already deleted the “global”, i realised that i dont need it.

Thanks, i will do!

11

It’s not just the global annotation, though. The entire computation appears to occur in global scope, which makes everything global.

12

I’m afraid the code is not working: c_dDG2 = eval(parsed_dDG2) will error since DG3 etc. have not been defined yet. Moving this line to after the declaration/update in the loop also doesn’t work, as DG3 is defined only in the scope of the loop, while eval evaluates in global scope.

In my opinion it would be better to provide your original working code (where you did eval in the loop if I understand correctly) and then optimise it for better performance, rather than to provide a non-working attempt at faster performance. In any case, don’t worry, I sufficiently understand what’s going on in the code.


By the way, your current MWE is also not exactly minimal: you could work with a single equation and replace all of the physics with a simple dummy calculation. I guess a more extensive code snippet could be useful for more general performance improvements (e.g. the global scope issues mentioned by @DNF ), but let’s for now stick to evaling your updating equations only once :).

Concretely, you could use

using Underscores  # for nicer piping syntax
(...)
eq_motion1 = "dDG2=2*l1*m1*cos(DG3)*cos(DG5)*DG1+DG4*DG6+Mg1, dDG4=l1*m1*cos(DG3)*tan(DG5)*DG4-l1*m1*cos(DG3)*cos(DG5)*DG6+Mg1, dDG6=2*l1*m1*cos(DG3)*sin(DG5)*DG4+DG2/DG6+Mg1"
eq_motion = @_ eq_motion1 |> split(__, ',') .|> split(__, '=') .|> __[2] |> map("(DG1, DG2, DG3, DG4, DG5, DG6, l1, m1, Mg1) -> " * _, __)
# e.g. eq_motion[1] == "(DG1, DG2, DG3, DG4, DG5, DG6, l1, m1, Mg1) -> 2*l1*m1*cos(DG3)*cos(DG5)*DG1+DG4*DG6+Mg1"
fun_dDG2, fun_dDG4, fun_dDG6 = eval.(Meta.parse.(eq_motion))

(...)
        # Calculate accelerations: update using new variable values
        dDG2 = fun_dDG2(DG1, DG2, DG3, DG4, DG5, DG6, l1, m1, Mg1)
        dDG4 = fun_dDG4(DG1, DG2, DG3, DG4, DG5, DG6, l1, m1, Mg1)
        dDG6 = fun_dDG6(DG1, DG2, DG3, DG4, DG5, DG6, l1, m1, Mg1)
        # (We could also collect these variables into a Tuple and splat it for more compact code)

(This does error in the second iteration, as in the third equation we divide by DG6 which is initialised as 0. Your full equations presumably won’t have this problem.)

13

Hi, thanks for the suggestion, the simulation now runs very good!

I currently have 1600 timesteps and the simulation finishes in ~1.5s, thats more than enough for me…

I also put the loop and a lot of other parts of my code into functions, that was a good idea.

2 Likes