FREQUENCY DOMAIN

[marcomanelli2]

BODE DIAGRAMS OF A SYSTEM OF n MASSES

Generic M, C and K case

using ControlSystems
using Plots
using ControlSystemsBase
using LinearAlgebra
using SparseArraysn=size(M,1)

Define the system matrices

A = spzeros(2n, 2n)
B = spzeros(2n, n)
CC = spzeros(n, 2n)

#First quarter stays as it is

#Second quarter
A[1:n,n+1:2n]=Matrix(1I,n,n)

#Third quarter of the A matrix (Stiffness)
A[n+1:2n,1:n]=-M\K

#Fourth quarter of the A matrix (Damping)
A[n+1:2n,n+1:2n]=-M\R

#B matrix
B[n+1:2n,1:n]=inv(M)

#C matrix
CC[1:n, 1:n]=Matrix(1I,n,n)

Transfer Function and Bode Diagrams

G = ss(A,B,CC,0)
fig=bodeplot(G)

[marcomanelli2]

EXAMPLE: n masses connected in series

using ControlSystems
using Plots
using ControlSystemsBase
using LinearAlgebra
using SparseArrays

Define the system parameters

m = [2.0, 2.0]   # masses (kg)
c = [0.2, 0.2, 0.2]  # damping coefficients (Ns/m)
k = [2.0, 2.0, 2.0]   # spring constants (N/m)
n = length(m)    # number of masses

Define the system matrices

A = spzeros(2n, 2n)
B = spzeros(2n, n)
CC = spzeros(n, 2n)
M=Diagonal(m);
C=Tridiagonal(-c[2:n], c[1:n]+c[2:n+1], -c[2:n]);
K=Tridiagonal(-k[2:n], k[1:n]+k[2:n+1], -k[2:n]);

#First quarter stays as it is

#Second quarter
A[1:n,n+1:2n]=sparse(1I,n,n)

#Third quarter of the A matrix (Stiffness)
A[n+1:2n,1:n]=Tridiagonal(k[2:n]./m[2:n], -(k[1:n]+k[2:n+1])./m, k[2:n]./m[1:n-1])

#Fourth quarter of the A matrix (Damping)
A[n+1:2n,n+1:2n]=Tridiagonal(c[2:n]./m[2:n], -(c[1:n]+c[2:n+1])./m, c[2:n]./m[1:n-1])

#B matrix
B[n+1:2n,1:n]=Diagonal(1 ./ m)

#C matrix
CC[1:n, 1:n]=Matrix(1I,n,n)

Transfer Function and Bode Diagrams

G = ss(A,B,CC,0)
fig=bodeplot(G)

Bode Diagrams.pdf

[marcomanelli2]

Finding the resonance peaks of a generic system

N-DIMENSION

insert N:

#e.g.:

N=3

using LinearAlgebra
using SparseArrays
using Plots
using LinearMaps
using IterativeSolvers
using LaTeXStrings

#e.g.: random values

m = rand(0.01:0.1:5,N)
c = rand(0.2:0.01:0.5,N+1)
k = rand(1:0.1:30,N+1)

Constructing the matrices

#e.g.: masses connected in series

M = Diagonal(m)
C = Tridiagonal(-c[2:N], c[1:N]+c[2:N+1], -c[2:N])
K = Tridiagonal(-k[2:N], k[1:N]+k[2:N+1], -k[2:N])

Studying system amplification through eigenvalues of the transfer function

(The biggest eigenvalue has been estimated through applying the inverse power method to find the smallest eigenvalue of the inverse transfer function)

w = LinRange(0, 20, 600)  #frequency range
lamb = spzeros(length(w))
eigvectors = rand(ComplexF64,length(w),N)

for j in 1:length(w)
    om = w[j]
    G = -om^(2)*M + im*C*om + K  #creating the inverse transfer function related to this frequency
    F = lu(G)
    Fmap = LinearMap{ComplexF64}((y, x) -> ldiv!(y, F, x), N, ismutating = true)
    lambda, x = invpowm(Fmap, shift = 0, tol = 1e-4, maxiter = 200)  #finding minimum eig in absolute value
    lambda = abs(lambda);
    eigvectors[j,:] = x
    lambda = 1/lambda  #finding eig of the transfer function
    lamb[j] = lambda  #saving the value
end

plot(w,lamb,xlabel=L"\omega",ylabel="System Amplification",label=L"$\lambda_{max}$", linewidth=1.5)