pid2df_sim.m

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

function [d_through_p]=pid2df_sim(simulate,pnummt,pdenmt,pdeadmt,mnummt,mdenmt,mdeadmt,...
   pdnummt,pddenmt,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 Pre_fil_den Pre_fil_num

 %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 col > 1 % multi-transfer functions
   
   if length(Epsilon) == 1
 		Epsilon(2)=Epsilon(1);
	end
   
   pdeadx=zeros(col,1); 
   mdeadx=zeros(col,1); 
   mnumT={[] []};
   mdenT={[] []};
   pnumT={[] []};
   pdenT={[] []};
   for i=1:col   
		%Model
		XnumMs{i}=mt2poly(mnummt{1,1,i},y);
		XdenMs{i}=mt2poly(mdenmt{1,1,i},y);
		[mdeadx(i),mnumT{i},mdenT{i}]=mt2dt(mdeadmt{1,1,i},y);
      %Process
      if isempty(pnummt{1,1}) | isempty(pdenmt{1,1})
         XnumPs{i}=XnumMs{i};
         XdenPs{i}=XdenMs{i};
         pnumT{i}=mnumT{i};
         pdenT{i}=mdenT{i};
         pdeadx(i)=mdeadx(i);
      else
         XnumPs{i}=mt2poly(pnummt{1,1,i},x);
			XdenPs{i}=mt2poly(pdenmt{1,1,i},x);
			[pdeadx(i),pnumT{i},pdenT{i}]=mt2dt(pdeadmt{1,1,i},x);
      end
           
   end %for i=1:col
   
   %Combind terms (normally two) to form approximated model with no dead time.
   disp('Approximated Model = tf(numMs,denMs) : in')
   if (length(XdenMs{1})==length(XdenMs{2}))& all(XdenMs{1}==XdenMs{2})
      temp1=conv(XnumMs{1},mnumT{1});
      temp1=tf(temp1,mdenT{1});
      temp2=conv(XnumMs{2},mnumT{2});
      temp2=tf(temp2,mdenT{2});
      temp3=tf(1,XdenMs{1});
      temp1=(temp1+temp2)*temp3;
   else
   	temp1=tf(XnumMs{1},XdenMs{1});
   	temp2=tf(mnumT{1},mdenT{1});
   	temp1=temp1*temp2;
   	temp2=tf(XnumMs{2},XdenMs{2});
   	temp3=tf(mnumT{2},mdenT{2});
   	temp2=temp2*temp3;
   	temp1=temp1+temp2;
   end
   [numMs,denMs]=tfdata(temp1,'v');
   tcf(tf(numMs,denMs))
   clear temp1 temp2 temp3
   
   %IMC 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 (with no dead time)
	[mq]=mq_gen(mnummt,mdenmt,mdeadmt,cfnummt,cfdenmt,orderf,y);
	[Csnum,Csden]=pid_tf2(y,numMs,denMs,1,1,denCs2,numCs2,...
   	denFs,denFs2,Tcanc,Epsilon,m,qf,inf);
   
   %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,m,qf,inf);
	% Disturbance lag------------------
	numPd=mt2poly(pdnummt,x);
	denPd=mt2poly(pddenmt,x);
   if any(real(roots(XdenPs{1}))>0) | any(real(roots(XdenMs{1}))>0) | any(real(roots(denPd))>0)
      p_unstable=1;
      if d_through_p==0
         d_through_p=1;
         warndlg('The plant is unstable --> Disturbance pass through the plant');
      end
   else
      p_unstable=0;
   end

else % if col = 1
   if length(Epsilon) == 1
      Epsilon(2)=Epsilon(1);
   end
   %Model
   [numMs]=mt2poly(mnummt{1,1},y);
   [denMs]=mt2poly(mdenmt{1,1},y);
   [mdeadx,mnumT,mdenT]=mt2dt(mdeadmt{1,1},y);
   Mdelay=mdeadx;
   %Process
   if isempty(pnummt{1,1}) | isempty(pdenmt{1,1})
      numPs=numMs;
      denPs=denMs;
      pdeadx=mdeadx;
      pnumT=mnumT;
      pdenT=mdenT;
   else
      [numPs]=mt2poly(pnummt{1,1},x);
      [denPs]=mt2poly(pdenmt{1,1},x);
      [pdeadx,pnumT,pdenT]=mt2dt(pdeadmt{1,1},x);
   end
   Pdelay=pdeadx;
   %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(mnummt{1},mdenmt{1},mdeadmt{1},cfnummt,cfdenmt,orderf,y);
   [Csnum,Csden]=pid_tf2(y,numMs,denMs,mnumT,mdenT,denCs2,numCs2,...
      denFs,denFs2,Tcanc,Epsilon,mq,orderf);
   % Controller part (for IMC simulation) & PID parameter
   pid2df(y,numMs,denMs,mnumT,mdenT,Csnum,Csden,denCs2,numCs2,...
      denFs,denFs2,Tcanc,Epsilon,order,mq,orderf);
	% Disturbance lag------------------
	numPd=mt2poly(pdnummt,x);
	denPd=mt2poly(pddenmt,x);
   if any(real(roots(denPs))>0) | any(real(roots(denMs))>0) | any(real(roots(denPd))>0)
      p_unstable=1;
      if d_through_p==0
         d_through_p=1;
         warndlg('The plant is unstable --> Disturbance pass through the plant');
      end
   else
      p_unstable=0;
   end
end
%-------------------------------------------




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

% Pre-filter part 
if strcmp(qd,'1')
   Qdnum=[1];
   Qdden=[1];
else
   [qd,Qdnum,Qdden]=qd_mat(Tcanc(2,:),mq,Epsilon(2),orderf,y);
end

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('Pre-Fillter in ');
tcf(tf(Pre_fil_num,Pre_fil_den));
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=15; 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;
   	if col > 1
      	if p_unstable
      		[t]=sim('imc2dfu2',time,options);
       	   temp1='* IMC (feedback form)'; 
   		else
      		[t]=sim('imc2df2',time,options);
       	   temp1='* IMC'; 
   		end
   	else
      	if p_unstable
      		[t]=sim('imc2dfu',time,options);
       	   temp1='* IMC (feedback form)'; 
   		else
      		[t]=sim('imc2df',time,options);
       	   temp1='* IMC'; 
      	end
   	end
      if plot_effort
         gh1=subplot(2,1,1);plot(t,simout,'m');
			title(['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(['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(['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(['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;