Skip to content

Instantly share code, notes, and snippets.

@pepijndevos
Last active March 12, 2023 18:58
Show Gist options
  • Save pepijndevos/043b0e95905014baa0106c3c134818bb to your computer and use it in GitHub Desktop.
Save pepijndevos/043b0e95905014baa0106c3c134818bb to your computer and use it in GitHub Desktop.
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment