How to find root locus of a polynomial

I need to find the root locus of a polynomial.

m_s = 1;    % Mass of the spring
m_r = 68;   % Mass of the rotor
k_r = 5.45e5; % Spring constant of the rotor
c_s = 663;  % Damping coefficient of the spring
c_r = 900;   % Damping coefficient of the rotor
k_s = 15000; % Spring constant of the spring

% Frequency range in rad/s
omega = 0:100:60000;

% Initialize arrays to store roots
roots_real = zeros(length(omega), 5);
roots_imag = zeros(length(omega), 5);

% Initialize array to store crossing points
crossing_omega = zeros(1, 5);
crossing_from_real_to_imag = zeros(1, 5);

% Calculate and store roots for each angular frequency
for i = 1:length(omega)
    a4_omega = c_s./m_s + k_r./c_r - 1i.*omega(i);
    a3_omega = k_r./m_r + k_r./m_s + k_s./m_s + c_s./m_s*(k_r./c_r - 1i.*omega(i));
    a2_omega = (k_r.*c_s) ./ (m_r.*m_s) + (k_r.*k_s) ./ (m_s.*c_r) - (k_r./m_r + k_r./m_s + k_r./m_s).*1i.*omega(i);
    a1_omega = k_r ./m_r * (k_s./m_s - 1i.*omega(i).*c_s./m_s);
    a0_omega = -1i.*omega(i).*k_s.*k_r ./(m_s.*m_r);

    a = [1, a4_omega, a3_omega, a2_omega, a1_omega, a0_omega];

    % Calculate roots
    system_roots = roots(a);
% Store real and imaginary parts of the roots
    roots_real(i, :) = real(system_roots);
    roots_imag(i, :) = imag(system_roots);

I would like to plot the root locus for this system as ω varies. Could you please assist me in implementing this in MATLAB? Specifically, I am looking for guidance on expressing the characteristic equation in terms of s and ω, and how to use MATLAB functions to generate the root locus plot

Answer 1

[Sulaymon Eshkabilov]

Here is how to get the roots plotted.

m_s = 1;    % Mass of the spring
m_r = 68;   % Mass of the rotor
k_r = 5.45e5; % Spring constant of the rotor
c_s = 663;  % Damping coefficient of the spring
c_r = 900;   % Damping coefficient of the rotor
k_s = 15000; % Spring constant of the spring

% Frequency range in rad/s
omega = 0:10000;

% Initialize arrays to store roots
roots_real = zeros(length(omega), 5);
roots_imag = zeros(length(omega), 5);

% Initialize array to store crossing points
crossing_omega = zeros(1, 5);
crossing_from_real_to_imag = zeros(1, 5);

% Calculate and store roots for each angular frequency
for i = 1:length(omega)
    a4_omega = c_s./m_s + k_r./c_r - 1i.*omega(i);
    a3_omega = k_r./m_r + k_r./m_s + k_s./m_s + c_s./m_s*(k_r./c_r - 1i.*omega(i));
    a2_omega = (k_r.*c_s) ./ (m_r.*m_s) + (k_r.*k_s) ./ (m_s.*c_r) - (k_r./m_r + k_r./m_s + k_r./m_s).*1i.*omega(i);
    a1_omega = k_r ./m_r * (k_s./m_s - 1i.*omega(i).*c_s./m_s);
    a0_omega = -1i.*omega(i).*k_s.*k_r ./(m_s.*m_r);

    a = [1, a4_omega, a3_omega, a2_omega, a1_omega, a0_omega];

    % Calculate roots
    system_roots = roots(a);
% Store real and imaginary parts of the roots
    roots_real(i, :) = real(system_roots);
    roots_imag(i, :) = imag(system_roots);
end
%%
figure
plot(roots_real(:,1), roots_imag(:,1), 'ro')
hold all
plot(roots_real(:,2), roots_imag(:,2), 'bo')
plot(roots_real(:,3), roots_imag(:,3), 'mx')
plot(roots_real(:,4), roots_imag(:,4), 'k*')
plot(roots_real(:,5), roots_imag(:,5), 'kd')
ylim([-20 50])
xlim([-16 6])
ylim([-1500 1e4])
xlim([-1050 200])
xlabel('Im(s)')
ylabel('Re(s)')

figure
plot(roots_real(:,1), roots_imag(:,1), 'ro')
hold all
plot(roots_real(:,2), roots_imag(:,2), 'bo')
plot(roots_real(:,3), roots_imag(:,3), 'mx')
plot(roots_real(:,4), roots_imag(:,4), 'k*')
plot(roots_real(:,5), roots_imag(:,5), 'kd')

ylim([-20 50])
xlim([-16 6])
xlabel('Im(s)')
ylabel('Re(s)')

[Supreeth D K]

Tq so much.

Answer 2

[Sulaymon Eshkabilov]

Step 1. Derive the transfer function of the system. E.g.: T = tf(1, [1 2 3])

Step 2. Use rlocus() or rlocusplot(). E.g.: rlocus(T)

m_s = 1;    % Mass of the spring
m_r = 68;   % Mass of the rotor
k_r = 5.45e5; % Spring constant of the rotor
c_s = 663;  % Damping coefficient of the spring
c_r = 900;   % Damping coefficient of the rotor
k_s = 15000; % Spring constant of the spring
% Frequency range in rad/s
omega = 0:100:6000;   % Smaller range
% Initialize arrays to store roots
roots_real = zeros(length(omega), 5);
roots_imag = zeros(length(omega), 5);
% Initialize array to store crossing points
crossing_omega = zeros(1, 5);
crossing_from_real_to_imag = zeros(1, 5);
% Calculate and store roots for each angular frequency
for i = 1:length(omega)
    a4_omega = c_s/m_s + k_r/c_r - omega(i);
    a3_omega = k_r/m_r + k_r/m_s + k_s/m_s + c_s/m_s*(k_r/c_r - omega(i));
    a2_omega = (k_r*c_s) / (m_r*m_s) + (k_r*k_s)/ (m_s*c_r) - (k_r/m_r + k_r/m_s + k_r/m_s)*omega(i);
    a1_omega = k_r /m_r * (k_s./m_s - omega(i)*c_s/m_s);
    a0_omega = -omega(i)*k_s*k_r/(m_s*m_r);
    a = [1, a4_omega, a3_omega, a2_omega, a1_omega, a0_omega];
    SYS = tf(1, a);
    rlocus(SYS), hold all
% % Calculate roots
%     system_roots = roots(a);
% % Store real and imaginary parts of the roots
%     roots_real(i, :) = real(system_roots);
%     roots_imag(i, :) = imag(system_roots);
end

[Supreeth D K]

But this is the plot they have obtained.