flocus.m

% This program plots the root loci of a fractional order system with
% transfer function
%
%             m/v       (m-1)/v
%         b  s   + b   s       + ... + b 
%          m        m-1                 0
% G(s) = --------------------------------- ,         (n>m),  
%             n/v       (n-1)/v
%         a  s   + a   s       + ... + a
%          n        n-1                 0
%
% on the first Riemann sheet for 0 < k < k_max. The vectors "num" and "den"
% are defined as follows:
% 
% num = [b   b    ... b ]
%         m   m-1      0
% den = [a   a    ... a ]
%         n   n-1      0
%
% The user must determine the values of "num", "den", "v", and "k_max". The
% value of "n" is limited to few hundred. 
%
% Example 1:
%
%   v = 2;
%   den = [1 -3 -2 2 12];
%   num = [1 -1];
%   k_max = 30;
%   
% Example 2:
%
%   v = 100;
%   den = [14994 zeros(1,33) 6009.5 zeros(1,96) 1.69];
%   num = [1];
%   k_max = 35;
%
% copyright (c) Farshad Merrikh-Bayat

clear all

v = 10;
den = [.8 zeros(1,22-9-1) .5 zeros(1,8) 1];
num = [1 zeros(1,11) 1];

k_max = 10;     % The root loci will be plotted for 0<k<k_max

%%%%%%%%%%%%%%%%%%%%%%%%%%%% main program %%%%%%%%%%%%%%%%%%%%%%%%%%
num = [zeros(1,length(den)-length(num)) num];

tol = 0.001;    % Smaller the value of "tol" more precise the plot. 
                % The value of tol = 0.001 is found to be suitable.

last_k = 0;
dk = .1;        % The initial step size for k. It may be automatically 
                % decreased if necessary. 
                
last_temp1 = roots(den);

k = last_k+dk;

figure(1)
clf
hold on

while k<k_max
    c = 1;
    k = last_k+dk;
    temp1 = roots(den+k*num);
    while (max(abs(temp1-last_temp1))>tol) & (c<2^20)
        c = 2*c;
        k = last_k+dk/c;
        temp1 = roots(den+k*num);
    end
    last_k = k;
    last_temp1 = temp1;
    sprintf('k=%f',k)
    s = 1;
    temp2 = [];
    for j=1:length(temp1)
        if angle(temp1(j))<=pi/v & angle(temp1(j))>-pi/v
            temp2(s) = temp1(j);
            s = s+1;
        end
    end
    if length(temp2)>0
        x = real(temp2.^v);
        y = imag(temp2.^v);
        h1 = plot(x,y,'k.');
        set(h1,'markersize',0.05)
    end
end  

temp3 = roots(den);
temp4 = [];
s = 1;
for k=1:length(temp3)
    if angle(temp3(k))<pi/v & angle(temp3(k))>-pi/v
        temp4(s) = temp3(k).^v;
        % The vector "temp4" contains the location of open-loop poles
        s = s+1;
    end
end
if length(temp4)>0
    h2 = plot(real(temp4),imag(temp4),'kx');
    set(h2,'markersize',10)
end

temp5 = roots(num);
temp6 = [];
s = 1;
for k=1:length(temp5)
    if angle(temp5(k))<pi/v & angle(temp5(k))>-pi/v
        temp6(s) = temp5(k).^v;
        % The vector "temp6" contains the location of open-loop zeros
        s = s+1;
    end
end

if length(temp5)>0 & length(temp6)>0 % To make sure that the system has at least one zero
    plot(real(temp6),imag(temp6),'ko');
end

grid
xlabel('Re\{s\}')
ylabel('Im\{s\}')
hold off