bicycle.jl

using ControlSystems
using Plots

# calculate the inertia from bifilar pendulum experiment
# https://www.youtube.com/watch?v=IdhV3lphRcc

L = 1.8 # length of the pendulum
D = 0.18 # distance between the wires
b = D/2
g = 9.81 # graviational constant
T = 2.01 # measured period
m = 0.468 # mass of the bicycle
Ix = m*g*b^2*T^2/(4*π^2*L)

L = 1.75 # length of the pendulum
D = 0.22 # distance between the wires
b = D/2
T = 1.38 # measured period
Iz = m*g*b^2*T^2/(4*π^2*L)

# A simple bike model
# https://lucris.lub.lu.se/ws/portalfiles/portal/4692388/625565.pdf

Vm = 2010/360 # measured velocity in rotations/s
cir = 0.08*π # wheel circumference
V = Vm*cir # velocity

a = 0.10 # distance from back wheel to CG
b = 0.18 # wheel base
h = 0.12 # height of CG
J = Ix

# For the inertia approximations you want h to represent some form of
# "square average" distance to the center of mass from the rotation axis.
# If you let h be the furthest distance from the axis, you'll over estimate the inertia
heff = sqrt(Ix/m)
aeff = sqrt(Iz/m)
D = m*aeff*heff

s = tf("s")
Pb = V*D/(b*J)*(s+m*V*h/D)/(s^2-m*g*h/J)
Ps = delay(0.01)
ζ = 0.8
ω = 17
# Pm = 1/s*tf([ω^2], [1, 2*ζ*ω, ω^2])
Pm = tf([1], [0.14, 1])

P = Ps * Pb * Pm
C = pid(10, 0.1, 0.08, Tf=0.001)
C2 = pid(6, 0.6, 0.06, Tf=0.001)
Cll = leadlinkat(45, 4)
nyquistplot([P*C, P*C2*Cll], Ms_circles=1.5, Mt_circles=1.5, ylims=(-2,5), xlims=(-5, 1), lab=["PID" "PID + leadlink"])

sys = feedback(P*C2*Cll)

plot(step(sys, 8), ylab="ϕ")

using JuliaSimControl, Plots
Ts = 0.01 # Discretization time
Tf = 25  # Simulation time

# Robustness constraints
Ms = 10   # Maximum allowed sensitivity function magnitude
Mt = Ms    # Maximum allowed complementary sensitivity function magnitude
Mks = 1000.0 # Maximum allowed magnitude of transfer function from process output to control signal, sometimes referred to as noise sensitivity.
w = 2π .* exp10.(LinRange(-2, 2, 200)) # frequency grid

prob = AutoTuningProblem(; P, Ms, Mt, Mks, w, Ts, Tf, metric = :IAE)

p0 = Float64[3, 0.3, 0.4, 0.001] # Initial parameter guess can be optionally supplied, kp, ki, kd, T_filter
res = solve(prob, p0)
plot(res)