pid.m

内模控制器（IMC）工具箱。包括参数整定、PID控制器参数转换等

function pid(numPs,denPs,numT,denT,numC,denC,numPinvs,denPinvs)

global K TauI TauD keps a1 a2 b1 b2 Alpha c1 c2 Kalf PID1lagden
global pid_text PID2lagden PIDrealden

% Calculates PID parameters from IMC design


% PID controller from numC/denC
%----------------------------------------------------------


REVnumC = fliplr(numC);			% reverse order of 'numC'
REVdenC = fliplr(denC);			% reverse order of 'denC'
REVnumPinvs = fliplr(numPinvs);		% reverse order of 'numPinvs'
REVdenPinvs = fliplr(denPinvs);		% reverse order of 'denPinvs'

keps=REVnumC(1)/REVdenC(2);
Kp=REVnumPinvs(1)/REVdenPinvs(1);

REVnumC = REVnumC/REVnumC(1);		% division by coeff. of linear term
REVdenC = REVdenC/REVdenC(2);   	% division by coeff. of linear term

if length(REVnumC) >= 3
 a1=REVnumC(2);
 a2=REVnumC(3);
elseif length(REVnumC) == 2
 a1=REVnumC(2);
 a2=0;
elseif length(REVnumC) <= 1
 a1=0;
 a2=0;
end

if length(REVdenC) >= 4
 b1=REVdenC(3);
 b2=REVdenC(4);
elseif length(REVnumC) == 3
 b1=REVdenC(3);
 b2=0;
elseif length(REVnumC) <= 2
 b1=0;
 b2=0;
end

pid_text='';
pid_text1='';


if  a1 ~= b1
 K    = keps*(a1-b1);
 TauI = a1-b1;
 TauD = (a2+b1^2-b2-a1*b1)/(a1-b1);
else
 pid_text='I cannot get a PID controller. Check entries.';
 say(pid_text);
 %save pid_cont pid_text
 return
end


% ideal PID controller (TauI, TauD positive)
%----------------------------------------------------------

numPIDi=keps*[TauI*TauD TauI 1];
denPIDi=[1 0];

disp(' ');
disp(' The ideal PID controller is:           ');
disp('                                        ');
disp('                    1                   ');
disp(' PI(s) = K*( 1 +  ------ + TauD*s)      ');
disp('                  TauI*s                ');
disp('				              ');
disp([' K    (controller gain)          = ', num2str(K)]);                        
disp([' TauI (integral time constant)   = ', num2str(TauI)]);
disp([' TauD (derivative time constant) = ', num2str(TauD)]);
disp(' ');

pid_text='';

pid_text=strvcat(pid_text, ' The ideal PID controller is: ');
pid_text=strvcat(pid_text, '                      1 ');
pid_text=strvcat(pid_text, ' PI(s) = K * (1 + ---------- + TauD*s ) ');
pid_text=strvcat(pid_text, '                    TauI*s ');
pid_text=strvcat(pid_text, 'where ');
pid_text0='';
pid_text0=['K                               = ', num2str(K)];
pid_text=strvcat(pid_text, pid_text0 );
pid_text0='';
pid_text0=['TauI                           = ', num2str(TauI)];
pid_text=strvcat(pid_text,pid_text0);
pid_text0='';
pid_text0=['TauD                          = ', num2str(TauD)];
pid_text=strvcat(pid_text, pid_text0);
pid_text=strvcat(pid_text, ' ********************************************************************************** ');
pid_text=strvcat(pid_text, '');



% lag PID controller (by adding p(s) )
%-------------------------------------

temp=[];
for i=1:length(denC)-1
 temp(i)=denC(i);
end
denC=temp;

[n1,d1]=polyder(numC,denC);
[n2,d2]=polyder(n1,d1);
[n3,d3]=polyder(n2,d2);

Alpha=-(n3(length(n3))/d3(length(d3))) / (n2(length(n2))/d2(length(d2)))/3;

w0=numC(length(numC))/denC(length(denC));

nw1=polyadd(conv(conv(n1,[Alpha 1]),denC),conv(Alpha*numC,d1));
dw1=conv(d1,denC);

w1=nw1(length(nw1))/dw1(length(dw1));

nw2=polyadd(conv(conv(n2,[Alpha 1]),d1),conv(2*Alpha*n1,d2));
dw2=conv(d1,d2);

w2=nw2(length(nw2))/dw2(length(dw2))/2;

Kalf=w0;
c1=w1/Kalf;
c2=w2/Kalf;

numPIDL1=Kalf*[c2 c1 1];
denPIDL1=conv([1 0],[Alpha 1]);
temp1=conv(denPs,denPIDL1);
temp2=conv(numPs,numPIDL1);
PID1lagden=polyadd(conv(temp1,denT),conv(temp2,numT));
temp1=roots(PID1lagden);
if any(real(temp1) > 0)
disp(' ');
disp(' The 1st order lag PID controller (blue) is unstable.');
disp(' ');
else
disp(' ');
disp(' The 1st order lag PID controller (blue) is:		');
disp('                             			');
disp('               k   (a2*s^2+a1*s+1)			');
disp(' PIDL1(s) = --- -----------------          	');
disp('              s     (Alpha*s+1)			');
disp('							');
disp([' where k                             = ', num2str(Kalf)]);                        
disp(['       a2                            = ', num2str(c2)]);
disp(['       a1                            = ', num2str(c1)]);
disp(['       Alpha                         = ', num2str(Alpha)]);
disp(' ');
if c1 < 0 | c2 < 0
disp(' The 1st order lag PID controller (blue) is unrealizable.');
end
end

pid_text=strvcat(pid_text, ' The 1st order lag PID controller (blue) is: ');
pid_text=strvcat(pid_text,'                      k   a2*s^2+a1*s+1 ');
pid_text=strvcat(pid_text, ' PIDL1(s) =  --- -------------------------------  ');
pid_text=strvcat(pid_text, '                     s    ( Alpha*a+1 ) ');
pid_text=strvcat(pid_text, ' where ');
pid_text0='';
pid_text0=['K                               = ', num2str(Kalf)];
pid_text=strvcat(pid_text, pid_text0 );
pid_text0='';
pid_text0=['a2                           = ', num2str(c2)];
pid_text=strvcat(pid_text,pid_text0);
pid_text0='';
pid_text0=['a1                          = ', num2str(c1)];
pid_text=strvcat(pid_text, pid_text0);
pid_text0=['Alpha                          = ', num2str(Alpha)];
pid_text=strvcat(pid_text, pid_text0);
pid_text=strvcat(pid_text, ' ********************************************************************************** ');


% lag PID controller (by dropping higher order terms)
%----------------------------------------------------

numPIDL=keps*[a2 a1 1];
denPIDL=conv([1 0],[b2 b1 1]);
temp1=conv(denPs,denPIDL);
temp2=conv(numPs,numPIDL);
PID2lagden=polyadd(conv(temp1,denT),conv(temp2,numT));
temp1=roots(PID2lagden);
if any(real(temp1) > 0)
disp(' ');
disp(' The 2nd order lag PID controller (green) is unstable.');
disp(' ');
else
disp(' ');
disp(' The 2nd order lag PID controller (green) is:		');
disp('                             			');
disp('            k   (a2*s^2+a1*s+1)			');
disp(' PIDL(s) = --- -----------------          	');
disp('            s   (b2*s^2+b1*s+1)			');
disp('							');
disp([' where k                             = ', num2str(keps)]);                        
disp(['       a2                            = ', num2str(a2)]);
disp(['       a1                            = ', num2str(a1)]);
disp(['       b2                            = ', num2str(b2)]);
disp(['       b1                            = ', num2str(b1)]);
disp(' ');
end


pid_text=strvcat(pid_text, ' The 2nd order lag PID controller (green) is: ');
pid_text=strvcat(pid_text, '                     k   a2*s^2+a1*s+1 ');
pid_text=strvcat(pid_text, ' PIDL(s) =  --- -------------------------------  ');
pid_text=strvcat(pid_text, '                     s   b2*s^2+b1*s+1) ');
pid_text=strvcat(pid_text, ' where ');
pid_text0='';
pid_text0=['k                               = ', num2str(keps)];
pid_text=strvcat(pid_text, pid_text0 );
pid_text0='';
pid_text0=['a2                           = ', num2str(a2)];
pid_text=strvcat(pid_text,pid_text0);
pid_text0='';
pid_text0=['a1                          = ', num2str(a1)];
pid_text=strvcat(pid_text, pid_text0);
pid_text0='';
pid_text0=['b2                           = ', num2str(b2)];
pid_text=strvcat(pid_text,pid_text0);
pid_text0='';
pid_text0=['b1                          = ', num2str(b1)];
pid_text=strvcat(pid_text, pid_text0);
pid_text=strvcat(pid_text, ' ********************************************************************************** ');


% realizable PID controler (TauI, TauD positive)
%----------------------------------------------------------

numPIDa1=keps*[(TauI*TauD+TauI*.05*TauD) (TauI+.05*TauD) 1];
denPIDa1=[.05*TauD 1 0];
temp1=conv(denPs,denPIDa1);
temp2=conv(numPs,numPIDa1);
PIDrealden=polyadd(conv(temp1,denT),conv(temp2,numT));
temp1=roots(PIDrealden);
if any(real(temp1) > 0)
disp(' ');
disp(' The realizable PID controller (red) is unstable.');
disp(' ');
else
disp(' ');
disp(' The realizable PID controller (red) is:	');
disp('                             			');
disp('                    1           TauD*s            ');
disp(' PIDA1(s) = K*(1 + ------ +  --------------- )    ');
disp('                   TauI*s    0.05*TauD*s+1        ');
disp('							');
disp([' K     (controller gain)          = ', num2str(K)]);                        
disp([' TauI  (integral time constant)   = ', num2str(TauI)]);
disp([' TauD  (derivative time constant) = ', num2str(TauD)]);
disp(' ');
end


pid_text=strvcat(pid_text, ' The realizable PID controller (red) is: ');
pid_text=strvcat(pid_text, '                               1            TauD*s  ');
pid_text=strvcat(pid_text, ' PIDA1(s) = K*(1 + --------- +  ------------------------- ) ');
pid_text=strvcat(pid_text, '                               TauI*s    0.05*TauD*s+1  ');
pid_text=strvcat(pid_text, ' where ');
pid_text0='';
pid_text0=[' K         (controller gain)              = ', num2str(K)];
pid_text=strvcat(pid_text, pid_text0 );
pid_text0='';
pid_text0=[' TauI    (integral time constant)    = ', num2str(TauI)];
pid_text=strvcat(pid_text,pid_text0);
pid_text0='';
pid_text0=[' TauD  (derivative time constant) = ', num2str(TauD)];
pid_text=strvcat(pid_text, pid_text0);
pid_text=strvcat(pid_text, ' ********************************************************************************** ');


%save pid_cont pid_text