Code covered by the BSD License

# AutotunerPID Toolkit

### William Spinelli (view profile)

• 1 file
• 4.22222

18 Mar 2004 (Updated )

Interactive simulation of a PID with autotuning.

bodePIDcompare(num,den,tau)
```function bodePIDcompare(num,den,tau)
%BODEPIDCOMPARE Comparison of Bode Diagrams with different autotuning
%   methods.
%
%   The function is typically called through the GUI in AutotunerPID.mdl
%   and provides a graphical comparison of four different autotuning
%   methods (STEP/ZN(OL), STEP/KT, STEP/IMC, RELAY/ZN(CL)) in frequency
%   domain.
%   By default the parameter Ms is set to 2 in KT methods, while the
%   parameter lambda is set 1/5 of the settling time of the plant.
%
%
%   BODEPIDCOMPARE(NUM,DEN,TAU) performs the comparative analysis directly
%   considering the plant described by the transfer function
%              NUM(s)
%      P(s) = -------- * exp(-TAU*s)
%              DEN(s)
%
%
%   Author:    William Spinelli (wspinell@elet.polimi.it)
%   \$Revision: 1.0 \$  \$Date: 2004/02/27 12:00:00 \$

global PIDPARAMETERS

% Bode diagram of the real plant model
P = tf(num,den,'IODelay',tau);
[m,f,w] = bode(P);
w = logspace(log10(w(1)),log10(w(end)),512);
s = j*w;

% compute step response
[y,t] = step(P);
Ts = min(real(abs(pole(P))))/50;
[y,t] = step(P,0:Ts:t(end));
Tf = t(end);

% identify a FOPDT model
modelFOPDT = idareas(y,1,Ts);

% identify the point of the frequency response with omega_n = -pi
[Gm,Pm,Wcg,Wcp] = margin(P);
modelPI.A = 2/Gm;
modelPI.T = 2*pi/Wcg;

% evaluate the frequency response of the plant
P = squeeze(freqresp(P,w)).';

% Create new figure
Lfig = figure(...
'Name',               'Bode diagram - Open loop function',...
'NumberTitle',        'off',...
Ffig = figure(...
'Name',               'Bode diagram - Complementary sensitivity',...
'NumberTitle',        'off',...
Sfig = figure(...
'Name',               'Bode diagram - Sensitivity',...
'NumberTitle',        'off',...

% Bode diagrams with different autotuner
% Ziegler & Nichols with Step identification
[K,Ti,Td,N,b] = pid_tuning(modelFOPDT,'ZN (OL)',[],'PID');
Rff = K*(b+1./(Ti*s));
Rfb = K*(1+1./(Ti*s)+(Td*s)./(1+Td*s/N));
L = Rfb.*P;
F = (Rff.*P)./(ones(1,length(P))+(Rfb.*P));
S = 1./(ones(1,length(P))+(Rfb.*P));
figure(Lfig)
subplot(211), semilogx(w,20*log10(abs(L)),'b')
subplot(212), semilogx(w,unwrap(angle(L))/pi*180,'b')
figure(Ffig)
subplot(211), semilogx(w,20*log10(abs(F)),'b')
subplot(212), semilogx(w,unwrap(angle(F))/pi*180,'b')
figure(Sfig)
subplot(211), semilogx(w,20*log10(abs(S)),'b')
subplot(212), semilogx(w,unwrap(angle(S))/pi*180,'b')

% Kappa-Tau with Ms = 2
% in order to obtain better performance
[K,Ti,Td,N,b] = pid_tuning(modelFOPDT,'KT',2,'PID');
Rff = K*(b+1./(Ti*s));
Rfb = K*(1+1./(Ti*s)+(Td*s)./(1+Td*s/N));
L = Rfb.*P;
F = (Rff.*P)./(ones(1,length(P))+(Rfb.*P));
S = 1./(ones(1,length(P))+(Rfb.*P));
figure(Lfig)
subplot(211), hold on, semilogx(w,20*log10(abs(L)),'r')
subplot(212), hold on, semilogx(w,unwrap(angle(L))/pi*180,'r')
figure(Ffig)
subplot(211), hold on, semilogx(w,20*log10(abs(F)),'r')
subplot(212), hold on, semilogx(w,unwrap(angle(F))/pi*180,'r')
figure(Sfig)
subplot(211), hold on, semilogx(w,20*log10(abs(S)),'r')
subplot(212), hold on, semilogx(w,unwrap(angle(S))/pi*180,'r')

% Internal Model Control with lambda = Tf/5
% this means that the closed loop model is made five time faster than the
% original (open loop model)
[K,Ti,Td,N,b] = pid_tuning(modelFOPDT,'IMC',Tf/5,'PID');
Rff = K*(b+1./(Ti*s));
Rfb = K*(1+1./(Ti*s)+(Td*s)./(1+Td*s/N));
L = Rfb.*P;
F = (Rff.*P)./(ones(1,length(P))+(Rfb.*P));
S = 1./(ones(1,length(P))+(Rfb.*P));
figure(Lfig)
subplot(211), hold on, semilogx(w,20*log10(abs(L)),'g')
subplot(212), hold on, semilogx(w,unwrap(angle(L))/pi*180,'g')
figure(Ffig)
subplot(211), hold on, semilogx(w,20*log10(abs(F)),'g')
subplot(212), hold on, semilogx(w,unwrap(angle(F))/pi*180,'g')
figure(Sfig)
subplot(211), hold on, semilogx(w,20*log10(abs(S)),'g')
subplot(212), hold on, semilogx(w,unwrap(angle(S))/pi*180,'g')

% Ziegler & Nichols with Relay identification
[K,Ti,Td,N,b] = pid_tuning(modelPI,'ZN (CL)',[],'PID');
Rff = K*(b+1./(Ti*s));
Rfb = K*(1+1./(Ti*s)+(Td*s)./(1+Td*s/N));
L = Rfb.*P;
F = (Rff.*P)./(ones(1,length(P))+(Rfb.*P));
S = 1./(ones(1,length(P))+(Rfb.*P));

figure(Lfig)
subplot(211), semilogx(w,20*log10(abs(L)),'m')
set(gca,'Position', [0.1514 0.4858 0.7536 0.3916],...
'XColor',[0.4 0.4 0.4],'YColor',[0.4 0.4 0.4],...
'FontSize',8,...
'XTickLabel',[])
title('Bode Diagram - Open loop function',...
'Color',[0 0 0],'FontSize',8)
ylabel('Magnitude (dB)',...
'Color',[0 0 0],'FontSize',8)
legend('STEP + ZN(OL)','STEP + KT','STEP + IMC','RELAY + ZN(CL)',4)
subplot(212), semilogx(w,unwrap(angle(L))/pi*180,'m')
set(gca,'Position',[0.1514 0.1100 0.7536 0.3472],...
'XColor',[0.4 0.4 0.4],'YColor',[0.4 0.4 0.4],...
'FontSize',8)
ylabel('Phase (deg)',...
'Color',[0 0 0],'FontSize',8)
'Color',[0 0 0],'FontSize',8)

figure(Ffig)
subplot(211), semilogx(w,20*log10(abs(F)),'m')
set(gca,'Position', [0.1514 0.4858 0.7536 0.3916],...
'XColor',[0.4 0.4 0.4],'YColor',[0.4 0.4 0.4],...
'FontSize',8,...
'XTickLabel',[])
title('Bode Diagram - Complementary sensitivity function',...
'Color',[0 0 0],'FontSize',8)
ylabel('Magnitude (dB)',...
'Color',[0 0 0],'FontSize',8)
legend('STEP + ZN(OL)','STEP + KT','STEP + IMC','RELAY + ZN(CL)',4)
subplot(212), semilogx(w,unwrap(angle(F))/pi*180,'m')
set(gca,'Position',[0.1514 0.1100 0.7536 0.3472],...
'XColor',[0.4 0.4 0.4],'YColor',[0.4 0.4 0.4],...
'FontSize',8)
ylabel('Phase (deg)',...
'Color',[0 0 0],'FontSize',8)
'Color',[0 0 0],'FontSize',8)

figure(Sfig)
subplot(211), semilogx(w,20*log10(abs(S)),'m')
set(gca,'Position', [0.1514 0.4858 0.7536 0.3916],...
'XColor',[0.4 0.4 0.4],'YColor',[0.4 0.4 0.4],...
'FontSize',8,...
'XTickLabel',[])
title('Bode Diagram - Sensitivity function',...
'Color',[0 0 0],'FontSize',8)
ylabel('Magnitude (dB)',...
'Color',[0 0 0],'FontSize',8)
legend('STEP + ZN(OL)','STEP + KT','STEP + IMC','RELAY + ZN(CL)',4)
subplot(212), semilogx(w,unwrap(angle(S))/pi*180,'m')
set(gca,'Position',[0.1514 0.1100 0.7536 0.3472],...
'XColor',[0.4 0.4 0.4],'YColor',[0.4 0.4 0.4],...
'FontSize',8)
ylabel('Phase (deg)',...
'Color',[0 0 0],'FontSize',8)