There are two major bugs in your code.
- num and den must be of the same length, which is probably not in your case. So, den+num.*k will not work if length(den) is NOT EQUAL to length(num). You can avoid it by using newnum instead of num. See the code for the definition of newnum.
- You must use the matrix as in qs(ii,jj) where ii and jj must be integers. So, if k=0.1, qs(k,:) will result in error.
- roots command may lead to another bug. See the code below.
I suggets the following code that still contains Bug 3 above. (See comments)
%% Plot ROOT LOCI without using the control Systems Tool Box command rlocus()
% my_loci(num, den, kend, usrtitle) plots the Root Loci of the open loop transfer function
% G(s)H(s)=N(s)/D(s)
% num are the coefficients of N(s)
% den are the coefficients of D(s)
%
% Written by P Vivekananda Shanmuganathan.
% Email: Vivek.Shanmuganathan@gmail.com
%
% Oct 24, 2019
% ========================
% CAUTION:
%
% There is a possible error in algorithm. The roots command may possibly
% change the order of the roots as the value of K is changed in each loop.
% This will lead to jump in the loci, i.e. roots of one locus will be
% mistaken for another.:-)
% For example, you can observe this error if you try with num = [3 10] and
% den = 1 1 1 1 ]
% ========================
function locus = my_loci(num, den, kend, usrtitle)
%% Create the matrix locus to hold root loci data.
locus = zeros(kend, length(den)-1);
%% Make dimensions of numerator same as the denominator and fill with zeros.
newnum = [zeros(1,length(den)-length(num)) num];
% If num=[3 1] and den=[1 1 1 1], newnum will be newnum=[0 0 3 1].
% This will make obtaining the coefficients of the characteristic
% polynomial as den+K*num .
kdata=0:kend; % Range of values for the gain K
for K=kdata % K is the Gain
% den+K*newnum is the characteristic polynomial.
% Roots of the chatacteristic polynimal are the poles for given K.
locus(K+1,:)=roots(den+K*newnum)';
end
%% Poles and zeros of the open loop
if length(den>1)
op = roots(den); % Open loop poles
opRe = real(op); % Real part
opIm = imag(op); % Imaginary part
end
if length(num>1)
oz = roots(num); % Open loop zeros
ozRe = real(oz); % Real part
ozIm = imag(oz); % Imaginary part
end
% Mark the open-loop poles and zeros and hold the graph.
clf; % Clear figure
plot(opRe, opIm,'x','LineWidth', 2, 'MarkerSize',15); hold;
plot(ozRe, ozIm,'o','LineWidth', 2, 'MarkerSize',8);
%% Plot the root loci.
% Each element of the locus contains a complex number.
% Plot command plots the imag parts against the real parts
plot(locus); % Plot the root loci
sgrid; % Display the special grid Radial lines indicate zeta and circles indicate w_n
title(usrtitle);
xlabel('Real Axis: Attenuation \sigma=\zeta\omega_n (seconds^-^1)')
ylabel('Imaginary Axis: Natural Frequency j\omega_n (seconds^-^1)')
axis equal % Scale same for both the axes
hold off
end
