pid_cascade2.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 K1 TauI1 TauD1 pid_text pid_text2
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)
   if all(XdenMs{1,1}==XdenMs{1,2})
      temp1=conv(XnumMs{1,1},XmnumT{1,1});
      temp1=tf(temp1,XmdenT{1,1});
      temp2=conv(XnumMs{1,2},XmnumT{1,2});
      temp2=tf(temp2,XmdenT{1,2});
      temp3=tf(1,XdenMs{1,1});
      temp1=(temp1+temp2)*temp3
   else
   	temp1=tf(XnumMs{1,1},XdenMs{1,1});
   	temp2=tf(XmnumT{1,1},XmdenT{1,1});
   	temp1=temp1*temp2;
   	temp2=tf(XnumMs{1,2},XdenMs{1,2});
   	temp3=tf(XmnumT{1,2},XmdenT{1,2});
   	temp2=temp2*temp3;
   	temp1=temp1+temp2
   end
   [numMs1,denMs1]=tfdata(temp1,'v');
   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,inf);
   [Csnum1,Csden1]=pid_tf3(conv(numMs,numMs1),conv(denMs,denMs1),denCs,numCs,denFs);
   %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);
   disp('Outter loop PID Controller:');
   [K1,TauI1,TauD1,pid_text2,err]=pid3(conv(numMs,numMs1),conv(denMs,denMs1),...
      Csnum1,Csden1,denCs,numCs);
   if err
      return
   end
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);
   numMs=conv(XnumMs{2,1},XmnumT{2,1});
   denMs=conv(XdenMs{2,1},XmdenT{2,1});
   numMs1=conv(XnumMs{1,1},XmnumT{1,1});
   denMs1=conv(XdenMs{1,1},XmdenT{1,1});
   [Csnum1,Csden1]=pid_tf3(conv(numMs,numMs1),conv(denMs,denMs1),...
      denCs,numCs,denFs);

	% 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);
   disp('Outter loop PID Controller:');
   [K1,TauI1,TauD1,pid_text2,err]=pid3(conv(numMs,numMs1),conv(denMs,denMs1),...
      Csnum1,Csden1,denCs,numCs);
   if err
      return
   end
end

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

denCs=conv(denCs,denFs);
denCs2=conv(denCs2,denFs2);
[qd,Qdnum,Qdden]=qd_mat(Tcanc(2,:),mq,Epsilon(2),orderf,y);
pid_text2=strvcat(pid_text2,' ');   
pid_text2=strvcat(pid_text2,'The outter loop PID controller must be cascade ');   
pid_text2=strvcat(pid_text2,'with the following lead-lag implementation (qqd invert):');   
pid_text2=strvcat(pid_text2,' ');   
pid_text2=strvcat(pid_text2,mtx2str(conv(denCs2,Qdden),2));   
pid_text2=strvcat(pid_text2,'-------------------------------------');   
pid_text2=strvcat(pid_text2,mtx2str(conv(numCs2,Qdnum),2));   
pid_text2=strvcat(pid_text2,' ');   
pid_text2=strvcat(pid_text2,'**********************************************************');   
pid_text=strvcat(pid_text,pid_text2);

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 (both inner & outter loop)'; 
      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;