sanko.m

国外人编写的一个分数阶系统的控制工具箱

function [G, J, handle] = sanko (w, gain, phase, Q, n, m, toler, NMaxIters)

% [G, J, handle] = sanko (w, gain, phase, Q, n, m, toler, NMaxIters)
% This function returns a model for a plant
% with a known frequency behaviour.
% The model is
% b(1)*s^(m*Q) + b(2)*s^((m-1)Q) + ... + b(end-2)*s^(2*Q) + b(end-1)*s^Q + b(end)
% -------------------------------------------------------------------------------
% a(1)*s^(n*Q) + a(2)*s^((n-1)Q) + ... + a(end-2)*s^(2*Q) + a(end-1)*s^Q + a(end)
% The function receives the gain in dB and the phase in degrees
% for frequencies, in rad/s, specified in w.
% Parameters Q, n and m must be supplied.
% Parameter toler is optional; two consecutive iterations resulting in
% vectors of parameters with norms closer than toler will stop the iterations.
% Parameter NMaxIters is optional and is the maximum number of iterations.
% A structure G is returned with fields num and den,
% containing vectors b and a respectively.
% J is the quadratic error per sampling frequency.
% If the last output argument exists, a handle for a Bode plot will be returned.
% Duarte Val閞io 2004

if nargin < 7 | isempty(toler), toler = 10^(-3-log(n+m)); end
if nargin < 8 | isempty(NMaxIters), NMaxIters = 100; end

if size(w, 1) < size(w, 2)
    w = w'; % now w is a column vector for sure
end
if size(gain, 1) < size(gain, 2)
    gain = gain'; % now w is a column vector for sure
end
if size(phase, 1) < size(phase, 2)
    phase = phase'; % now w is a column vector for sure
end

R = real(10.^(gain/20) .* exp(j * deg2rad(phase)));
I = imag(10.^(gain/20) .* exp(j * deg2rad(phase)));
weight = 1;
coefficientsOld = zeros(n+m+1, 1);
iters = 0; % number of iterations performed
J = [];

while 1
    
    A = [];
    for l = 0:m
        for c = 0:m
            A(l+1,c+1) = sum((-real((j*w).^(l*Q)).*real((j*w).^(c*Q))...
                -imag((j*w).^(l*Q)).*imag((j*w).^(c*Q))) .* weight);
        end
    end
    
    B = [];
    for l = 0:m
        for c = 1:n
            B(l+1,c) = sum((real((j*w).^(l*Q)).*real((j*w).^(c*Q)).*R...
                +imag((j*w).^(l*Q)).*real((j*w).^(c*Q)).*I...
                -real((j*w).^(l*Q)).*imag((j*w).^(c*Q)).*I...
                +imag((j*w).^(l*Q)).*imag((j*w).^(c*Q)).*R) .* weight);
        end
    end
    
    C = [];
    for l = 1:n
        for c = 0:m
            C(l,c+1) = sum(( (imag((j*w).^(l*Q)).*I - real((j*w).^(l*Q)).*R) .* real((j*w).^(c*Q))...
                +(-real((j*w).^(l*Q)).*I - imag((j*w).^(l*Q)).*R) .* imag((j*w).^(c*Q)) ) .*...
                weight);
        end
    end
    
    D = [];
    for l = 1:n
        for c = 1:n
            D(l,c) = sum((real((j*w).^(l*Q)).*real((j*w).^(c*Q)).*(R.^2+I.^2)...
                +imag((j*w).^(l*Q)).*imag((j*w).^(c*Q)).*(R.^2+I.^2)) .* weight);
        end
    end
    
    e = [];
    for l = 0:m
        e(l+1,1) = sum((-real((j*w).^(l*Q)).*R -imag((j*w).^(l*Q)).*I) .* weight);
    end
    
    f = [];
    for l = 1:n
        f(l,1) = sum((-real((j*w).^(l*Q)).*(R.^2+I.^2)) .* weight);
    end
    
    coefficients = [A B; C D] \ [e; f];
    if sum(isinf(coefficients)) | sum(isnan(coefficients))
        coefficients = coefficientsOld; % avoids inf and NaN
    end
    iters = iters+1;
    G.num = coefficients(m+1:-1:1);
    G.den = [coefficients(end:-1:m+2); 1];

    % index J
    tempNum = 0;
    for c = 1:length(G.num)
        tempNum = tempNum + G.num(c) * (j*w).^((m-c+1)*Q);
    end
    tempDen = 0;
    for c = 1:length(G.den)
        tempDen = tempDen + G.den(c) * (j*w).^((n-c+1)*Q);
    end
    J = [J; iters sum((abs(10.^(gain/20).*exp(j*deg2rad(phase)) - tempNum./tempDen)).^2) / length(w)];

    if (norm(coefficientsOld - coefficients) < toler) | (iters > NMaxIters)
        break
    else
        weight = 1 ./ (abs(tempDen).^2);
        coefficientsOld = coefficients;
    end
end

% Bode plots
if nargout > 2
    wNew = logspace(floor(log10(w(1))), ceil(log10(w(end))), max(50, length(w)));
    tempNum = 0;
    for c = 1:length(G.num)
        tempNum = tempNum + G.num(c) * (j*wNew).^((m-c+1)*Q);
    end
    tempDen = 0;
    for c = 1:length(G.den)
        tempDen = tempDen + G.den(c) * (j*wNew).^((n-c+1)*Q);
    end
    gainNew = 20 * log10(abs(tempNum ./ tempDen));
    phaseNew = rad2deg(angle(tempNum ./ tempDen));
    handle = figure;
    subplot(2,1,1)
    semilogx(w,gain,'.', wNew,gainNew)
    subplot(2,1,2)
    semilogx(w,phase,'.', wNew,phaseNew)
end