pid_cascade.m

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

function pid_cascade(simulate,Xpnummt,Xpdenmt,Xpdeadmt,Xmnummt,Xmdenmt,Xmdeadmt,...
   Xpdnummt,Xpddenmt,cnummt,cdenmt,cfnummt,cfdenmt,order,orderf,x,y,Epsilon,...
   subs,slbs,m,qf,qd,Tcanc,d_through_p,col)
% Setpoint and disturbance step change for a plant without saturation bounds
%------------------------------------------------------------------------
   
global K TauI TauD PID2lagden PID_FIden  keps a2 a1 b2 b1 PID_APXden keep gh1 gh2
global PI_lagden PI_deng PID_realden KI TauX TauY time st_size plot_setpoint plot_disturb
global plot_effort setpoint_Enable disturb_Enable effort_Enable plot_IMC PI_Enable
global plot_real_PID plot_lag_PID plot_FI plot_PID_APX plot_PI_lag plot_PI PI_lag_Enable
global IMC_Enable real_PID_Enable lag_PID_Enable FI_Enable PID_APX_Enable  old_handle
global XnumMs XdenMs mdeadx XnumPs XdenPs pdeadx numPd denPd Pdelay Mdelay numCs denCs
global numCs2 denCs2 Csnum Csden numMs denMs numPs denPs Setpoint Disturb Qdnum Qdden
global Xmdeadx Xpdeadx Pre_fil_den Pre_fil_num pid_text pid_text1
Xpnummt=shiftdim(Xpnummt,2)';
Xpdenmt=shiftdim(Xpdenmt,2)';
Xpdeadmt=shiftdim(Xpdeadmt,2)';
Xmnummt=shiftdim(Xmnummt,2)';
Xmdenmt=shiftdim(Xmdenmt,2)';
Xmdeadmt=shiftdim(Xmdeadmt,2)';

%Filter
if order==0
	denFs=1;
else
denF=[];
for i=1:order
 denF(i,1)=Epsilon(1);
 denF(i,2)=1;
end

denFs=denF(1,:);
for i=2:order
 denFs=conv(denFs,denF(i,:));
end
end
if Epsilon(1)==0
 denFs=1;
end

if orderf==0
	denFs2=1;
else
denF2=[];
for i=1:orderf
 denF2(i,1)=Epsilon(2);
 denF2(i,2)=1;
end
denFs2=denF2(1,:);
for i=2:orderf
 denFs2=conv(denFs2,denF2(i,:));
end
end
if Epsilon(2)==0
 denFs2=1;
end

   if isempty(Xmnummt{1,1}) | isempty(Xmdenmt{1,1}) | isempty(cfnummt) | ...
   	isempty(cfdenmt) 
 		errordlg('Not enough data entered.')
 		return
    end
    XnumMs={[0] [0];[0] [0]};
    XdenMs={[1] [1]; [1] [1]};
    XnumPs={[0] [0];[0] [0]};
    XdenPs={[1] [1]; [1] [1]};
    Xmdeadx={[eps] [eps];[eps] [eps]};
    Xpdeadx={[eps] [eps];[eps] [eps]};
    XmnumT={[1] [1]; [1] [1]};
    XmdenT={[1] [1]; [1] [1]};
    XpnumT={[1] [1]; [1] [1]};
    XpdenT={[1] [1]; [1] [1]};
   for i=1:2
      for j=1:col
			%Model
			XnumMs{i,j}=mt2poly(Xmnummt{i,j},y);
			XdenMs{i,j}=mt2poly(Xmdenmt{i,j},y);
			[Xmdeadx{i,j},XmnumT{i,j},XmdenT{i,j}]=mt2dt(Xmdeadmt{i,j},y);
         Xmdeadx{i,j}=max([Xmdeadx{i,j} eps]);
      	%Process
      	if isempty(Xpnummt{1,1}) | isempty(Xpdenmt{1,1})
         	XnumPs{i,j}=XnumMs{i,j};
         	XdenPs{i,j}=XdenMs{i,j};
         	XpnumT{i,j}=XmnumT{i,j};
         	XpdenT{i,j}=XmdenT{i,j};
         	Xpdeadx{i,j}=mdeadx{i,j};
      	else
         	XnumPs{i,j}=mt2poly(Xpnummt{i,j},x);
				XdenPs{i,j}=mt2poly(Xpdenmt{i,j},x);
				[Xpdeadx{i,j},XpnumT{i,j},XpdenT{i,j}]=mt2dt(Xpdeadmt{i,j},x);
	         Xpdeadx{i,j}=max([Xpdeadx{i,j} eps]);
      	end
           
   	end %for j=1:col
	end % for i
if col > 1
   %Combind terms of the inner loop model(normally two) 
   %to form an approximated model with no dead time.
   disp('Approximated Model = tf(numMs,denMs) :')
   if all(XdenMs{2,1}==XdenMs{2,2})
      temp1=conv(XnumMs{2,1},XmnumT{2,1});
      temp1=tf(temp1,XmdenT{2,1});
      temp2=conv(XnumMs{2,2},XmnumT{2,2});
      temp2=tf(temp2,XmdenT{2,2});
      temp3=tf(1,XdenMs{2,1});
      temp1=(temp1+temp2)*temp3
   else
   	temp1=tf(XnumMs{2,1},XdenMs{2,1});
   	temp2=tf(XmnumT{2,1},XmdenT{2,1});
   	temp1=temp1*temp2;
   	temp2=tf(XnumMs{2,2},XdenMs{2,2});
   	temp3=tf(XmnumT{2,2},XmdenT{2,2});
   	temp2=temp2*temp3;
   	temp1=temp1+temp2
   end
   [numMs,denMs]=tfdata(temp1,'v');
   tf(numMs,denMs)
   clear temp1 temp2 temp3
   
   %IMC outter loop Controller
	[denCs]=mt2poly(cnummt,y);
	[numCs]=mt2poly(cdenmt,y);
   
   % Inner loop feedback controller q(e) 
   [denCs2]=mt2poly(cfnummt,y);
	[numCs2]=mt2poly(cfdenmt,y);

   %Approximated feed-back form controller (with no dead time)
	[mq]=mq_gen(Xmnummt(2,:),Xmdenmt(2,:),Xmdeadmt(2,:),cfnummt,cfdenmt,orderf,y);
	[Csnum,Csden]=pid_tf2(y,numMs,denMs,1,1,denCs2,numCs2,...
   	denFs,denFs2,Tcanc,Epsilon,mq,orderf);
   
   %Compute PID from approximated feed-back form controller 
   %with no dead time (numT=1,denT=1)
	pid2df(y,numMs,denMs,1,1,Csnum,Csden,denCs2,numCs2,...
   	denFs,denFs2,Tcanc,Epsilon,order,mq,inf);
else
	%Controller
	[denCs]=mt2poly(cnummt,y);
	[numCs]=mt2poly(cdenmt,y);

	% Feedback path controller q(e) 
	[denCs2]=mt2poly(cfnummt,y);
	[numCs2]=mt2poly(cfdenmt,y);

	%Approximated feed-back form controller 
	[mq]=mq_gen(Xmnummt{2,1},Xmdenmt{2,1},Xmdeadmt{2,1},cfnummt,cfdenmt,orderf,y);
   [Csnum,Csden]=pid_tf2(y,XnumMs{2,1},XdenMs{2,1},XmnumT{2,1},XmdenT{2,1},...
      denCs2,numCs2,denFs2,denFs2,Tcanc,Epsilon,mq,orderf);

	% Controller part (for IMC simulation) & PID parameter
   pid2df(y,XnumMs{2,1},XdenMs{2,1},XmnumT{2,1},XmdenT{2,1},Csnum,Csden,...
      denCs2,numCs2,denFs,denFs2,Tcanc,Epsilon,order,mq,orderf);
end
   

%-------------------------------------------



denCs=conv(denCs,denFs);
denCs2=conv(denCs2,denFs2);

% Pre-filter part 
[qd,Qdnum,Qdden]=qd_mat(Tcanc(2,:),mq,Epsilon(2),orderf,y);

Pre_fil_num=conv(numCs,denCs2);
Pre_fil_num=conv(Pre_fil_num,Qdden);
Pre_fil_den=conv(denCs,numCs2);
Pre_fil_den=conv(Pre_fil_den,Qdnum);

[Pre_fil_num,Pre_fil_den]=z_p_cancel(Pre_fil_num,Pre_fil_den,1);
disp('qr(s)/qqd(s) in ');
[temp,temp1,temp2,temp3]=tcf(tf(Pre_fil_num,Pre_fil_den));
pid_text1='qr(s)/qqd(s) in ';
pid_text1=strvcat(pid_text1,temp3);
clear temp temp1 temp2 temp3
Setpoint=1;
Disturb=0;
st_size=1;
TauX=TauI+abs(TauD);
TauY=abs(TauD);

%-------------------------------------
if ~simulate
   return
end
options = simset('SrcWorkspace','current');
if isempty(keep) | ~keep
if exist('old_handle') & ~isempty(old_handle) & old_handle==gcf
   figure(old_handle);
   if ~isempty(x)
      set(old_handle,'Name',['Step responses for process with current parameter x = [ ',...
            num2str(x),' ]']);
   else
      set(old_handle,'Name','Step responses for process with perfect model');
   end
   clf
 uicontrol(old_handle,...
  'Style','push',...
  'Position',[550 170 85 25],...
  'String','Change Process',...
  'Callback',['change_P(''main'',x,Epsilon{1,1});']);
 uicontrol(old_handle,...
  'Style','push',...
  'Position',[555 140 75 25],...
  'String','Add&Remove',...
  'Callback',['add_remove2(''add'');']);
uicontrol(old_handle,...
  'Style','push',...
  'Position',[555 110 75 25],...
  'String','Change time',...
  'Callback',['add_remove2(''Time'');']);
uicontrol(old_handle,...
  'Style','push',...
  'Position',[565 80 55 25],...
  'String','Close',...
  'Callback','close');
uicontrol(old_handle,...
  'Style','push',...
  'Position',[565 50 55 25],...
  'String','Help',...
  'Callback','sayhelp(4270)');
else
 time=50; st_size=1; plot_setpoint=1; plot_disturb=1; plot_effort=0;
 setpoint_Enable='on'; disturb_Enable='on'; effort_Enable='off';
 
 plot_IMC=1; plot_real_PID=1; plot_lag_PID=0; plot_FI=0;
 plot_PID_APX=0; plot_PI_lag=0; plot_PI=0;
 
 IMC_Enable='on'; real_PID_Enable='on'; lag_PID_Enable='on'; FI_Enable='on'; 
 PID_APX_Enable='on'; PI_lag_Enable='on'; PI_Enable='on';
 if ~isempty(x)
    temp=['Step responses for process with current parameter x = [ ',...
          num2str(x),' ]'];
 else
    temp='Step responses for process with perfect model';
 end
 f=crfig(20,40,634,410,temp,'k','figure','on');
 		whitebg('white');
 uicontrol(f,...
  'Style','push',...
  'Position',[550 170 85 25],...
  'String','Change Process',...
  'Callback',['change_P(''main'',x,Epsilon{1,1});']);
 uicontrol(f,...
  'Style','push',...
  'Position',[555 140 75 25],...
  'String','Add&Remove',...
  'Callback',['add_remove2(''add'');']);
uicontrol(f,...
  'Style','push',...
  'Position',[555 110 75 25],...
  'String','Change time',...
  'Callback',['add_remove2(''Time'');']);
uicontrol(f,...
  'Style','push',...
  'Position',[565 80 55 25],...
  'String','Close',...
  'Callback','close');
uicontrol(f,...
  'Style','push',...
  'Position',[565 50 55 25],...
  'String','Help',...
  'Callback','sayhelp(4270)');
end
else % keep
   if isempty(old_handle) | ~ishandle(old_handle)
      errordlg('The original figure was gone');
      return
   end
   figure(old_handle)
   if ~isempty(x)
      set(old_handle,'Name',['Step responses for process with current parameter x = [ ',...
            num2str(x),' ]']);
   else
      set(old_handle,'Name','Step responses for process with perfect model');
   end
end % keep

Satub=str2num(subs);
Satlb=str2num(slbs);
Setpoint=1;
Disturb=0;
y_scale=-0.2;
y_position=0.2;
x_position=0.5*time;

temp=Epsilon;
axis('auto');
if (plot_IMC)
   if plot_setpoint
      Setpoint=1;
   	Disturb=0;
  		[t]=sim('imccascade',time,options);
  	   temp1='* IMC'; 
      if plot_effort
         gh1=subplot(2,1,1);plot(t,simout,'m');
			title(['Cascade system Setpoint Step Response when Epsilon = ',num2str(temp)])
   		xlabel('Time'); ylabel('Output'); 
   		h=text(x_position,y_position,temp1);
   		set(h,'color','m');
   		y_position=y_position-y_scale;
   		grid on;
   		hold on
         
         gh2=subplot(2,1,2);plot(t,effort,'m');
			title(['Cascade system Control efforts when Epsilon = ',num2str(temp)])
   		xlabel('Time'); ylabel('Output'); 
   		grid on;
         hold on
      elseif plot_disturb
      	gh1=subplot(2,1,1);plot(t,simout,'m');
			title(['Cascade system Setpoint Step Response when Epsilon = ',num2str(temp)])
   		xlabel('Time'); ylabel('Output'); 
   		h=text(x_position,y_position,temp1);
   		set(h,'color','m');
   		y_position=y_position-y_scale;
   		grid on;
         hold on
      else
         gh1=plot(t,simout,'m');
			title(['Cascade system Setpoint Step Response when Epsilon = ',num2str(temp)])
   		xlabel('Time'); ylabel('Output'); 
   		h=text(x_position,y_position,temp1);
   		set(h,'color','m');
   		y_position=y_position-y_scale;
   		grid on;
         hold on
      end
   end
   if plot_disturb
   	Setpoint=0;
   	Disturb=1;
  		[t]=sim('imccascade',time,options);
  	   temp1='* IMC'; 
      if plot_effort
      	gh1=subplot(2,1,1);plot(t,simout,'m');
			title(['Cascade system Disturbance Step Response when Epsilon = ',num2str(temp)])
   		xlabel('Time'); ylabel('Output');