Flappy bird model of path and speed

In this quick example, we walk through a simple set of constraints on a 2-DOF point,

$$y = sin\pi x$$
$$\dot{x} = 2$$
,

where the postion is (x, y) and its velocity is $(\dot{x},~\dot{y})$. The point will travel forward at a constant speed, but follow a sine-wave path.

using Symbolics, DifferentialEquations, ModelingToolkit, Plots

using ModelingToolkit: t_nounits as t, D_nounits as D

import Latexify

Setting up functions and constants

In this example, we want to define the constraints as a Symbolic function. That way, we can take analytical derivatives and Jacobians. Here, set up the @variables and the function, C(q, t)

Note: x, y, vx, vy are all functions of time, but because we are not solving a set of differential equations we want to keep them as indepenent variables e.g. @variables xnot@variables x(t).

begin
    @variables x, y, vx, vy
    q = [x, y]
    dq = [vx, vy]
    C(q, t) = [y - sin(pi*x);
              x - 2*t]
end

C (generic function with 1 method)

Here, we can see the two constraints on motion written as the constraint equation, $C(q,~t)$ and its Jacobian, $\frac{\partial C}{\partial q}$. Each of these printed equation arrays are formed with the Julia Computer Algebra System so they can be rearranged and used in algebraic+calculus expressions.

C(q, 0.1)

$$\begin{equation} \left[ \begin{array}{c} y - \sin\left( x ~ \pi \right) \\ -0.2 + x \\ \end{array} \right] \end{equation}$$

dCdq(q, t) = Symbolics.jacobian(C(q, t), q)

dCdq (generic function with 1 method)

dCdq(q, t)*D.(q) .~ -Symbolics.derivative(C(q,t), t)

$$\begin{align} \frac{\mathrm{d}}{\mathrm{d}t} ~ y - \cos\left( x ~ \pi \right) ~ \pi ~ \frac{\mathrm{d}}{\mathrm{d}t} ~ x &= 0 \\ \frac{\mathrm{d}}{\mathrm{d}t} ~ x &= 2 \end{align}$$

Solving for actual positions and velocities

The equations are nice to see and check the math works out, but the goal is to have real, numerical coordinates, x-y, and real velocities, $v_x,~v_x$.

I break this down into three steps:

1. set up variables

N is my number of timesteps
tN is my range of time, here from 0-4 s with N steps
q_sol and dq_sol are initialized vectors of position and velocity, respectively
partialCt is the partial time derivative, $\frac{\partial C}{\partial t}$, of $C(q,t)$

2. solve for positions and velocities

There are some hoops to jump through to go from CAS to numerical results. Julia has a robust system of solvers for this task if we can match the syntax.

The equations we are solving are set up in

eqs = 0 .~ C(q, t0)

where t0 is the current time to solve the constraint equations and q is still a set of CAS variables, [x, y]. Then, we create the @named constraint equations with

@named constraint_eqn = NonlinearSystem(eqs, q, [])

where NonlinearSystem is a function in ModelingToolkit that defines a system at a given point in time. The complete function signals that the system is ready to solve. The values of x and y are saved in sol. Then, the q_sol[i,:] row is given the values of sol.

The solution for the velocities is a linear algebra equation, $\frac{\partial C}{\partial t}\dot{q} = -C_t$

$$\frac{\mathrm{d}}{\mathrm{d}t} \cdot y - \cos\left( x \cdot \pi \right) \cdot \pi \cdot \frac{\mathrm{d}}{\mathrm{d}t} \cdot x = 0$$

$$\frac{\mathrm{d}}{\mathrm{d}t} \cdot x = 2$$

This is where the conversion from CAS to numerical starts. The numerical values of x and y are saved in a Dict and used in the Symbolics.substitute function,

lhs = Symbolics.substitute(dCdq(q, t), sol_dict, fold=Val(true))
rhs = Symbolics.substitute(-partialCt, sol_dict, fold=Val(true))

where the fold=Val(true) turns the exact CAS statement into an array with a distinct value. Then,

solv = Symbolics.value.(lhs\rhs)

solves the Linear algebra expression for $\dot{q}$ and saves it as solv. These values are stored in dq_sol[i, :].

begin
    N = 100
    tN = LinRange(0, 4, N)
    q_sol = zeros(length(tN), 2)
    dq_sol = zeros(length(tN), 2)
    partialCt = Symbolics.derivative.(C(q,t), t)
    
    local sol = [0, 0]	
    
    for (i, t0) in enumerate(tN)
        
        eqs =0 .~ C(q,t0)
        @named constraint_dae = NonlinearSystem(eqs, q, [])
        complete_dae = complete(constraint_dae)
        u0 = [x => sol[1], y => sol[2]]
        prob = NonlinearProblem(complete_dae, u0)
        sol = solve(prob)
        q_sol[i,:] .= sol
    
        sol_dict = Dict(x =>sol[q][1], y => sol[q][2],t => t0)
        lhs = Symbolics.substitute(dCdq(q, t), sol_dict, fold=Val(true))
        rhs = Symbolics.substitute(-partialCt, sol_dict, fold=Val(true))
        solv = Symbolics.value.(lhs\rhs)
        dq_sol[i,:] .= solv
    end
end

Wrapping up

The process here depends upon identifying constraints and degrees of freedom carefully, but it can scale to as many constraints and degrees of freedom as our system might need. The code will change a little, but if we can define a system where #DOFs == # constraints, there are no differential equations and the problem is a set of nonlinear + linear algebra equations at each step in time.

Plus, we now have a Flappybird velocity profiler ;)

plot_path_with_velocity! (generic function with 1 method)

Built with Julia 1.12.5 and

DifferentialEquations 7.17.0
Latexify 0.16.10
ModelingToolkit 11.10.0
Plots 1.41.6
Symbolics 7.17.0

To run this tutorial locally, download this file and open it with Pluto.jl.

Dynamics