pid1df_sim.m

Software for design and tuninig of SISO and MIMO contol systems

function pid1df_sim(simulate,pnummt,pdenmt,pdeadmt,mnummt,mdenmt,mdeadmt,pdnummt,pddenmt,cnummt,cdenmt,...
   order,x,y,Epsilon,subs,slbs,d_through_p,col)
% Setpoint and disturbance step change for a plant without saturation bounds
%------------------------------------------------------------------------

global K TauI TauD keps a1 a2 b1 b2 Alpha c1 c2 Kalf PID1lagden PID2lagden PIDrealden 
global time st_size plot_setpoint plot_disturb plot_effort setpoint_Enable disturb_Enable 
global effort_Enable plot_IMC plot_real_PID plot_1_lag_PID plot_2_lag_PID IMC_Enable 
global real_PID_Enable lag1_PID_Enable lag2_PID_Enable old_handle keep
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

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

denFs=denF(1,:);
for i=2:order
 denFs=conv(denFs,denF(i,:));
end
end
if Epsilon==0
 denFs=1;
end
if col > 1 % multi-transfer functions
   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);
      mdeadx(i)=mdeadx(i)+eps;
      %Process
      if isempty(pnummt{1,1,1}) | isempty(pdenmt{1,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);
         pdeadx(i)=pdeadx(i)+eps;
      end
           
   end %for i=1:col
   
   %Combind two terms to form approximated model with no dead time.
   disp('Approximated Model = tf(numMs,denMs) :')
   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');
   if any(roots(denMs)> 0)
      unstable=1;
   else
      unstable=0;
   end
   tcf(tf(numMs,denMs))
   clear temp1 temp2 temp3
   
   %IMC Controller
	[denCs]=mt2poly(cnummt,y);
	[numCs]=mt2poly(cdenmt,y);
   
   %Approximated feed-back form controller (with no dead time)
   [Csnum,Csden]=pid_tf(numMs,denMs,1,1,denCs,numCs,denFs);
   
   %Compute PID from approximated feed-back form controller 
   %with no dead time (numT=1,denT=1)
	pid(numMs,denMs,1,1,Csnum,Csden,denCs,numCs);

else % if col = 1
   %Model
   [numMs]=mt2poly(mnummt{1,1},y);
   [denMs]=mt2poly(mdenmt{1,1},y);
   if any(roots(denMs)> 0)
      unstable=1;
   else
      unstable=0;
   end
   [mdeadx,mnumT,mdenT]=mt2dt(mdeadmt{1,1},y);
   %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
   %IMC Controller
   [denCs]=mt2poly(cnummt,y);
   [numCs]=mt2poly(cdenmt,y);
   %Approximated feed-back form controller
   [Csnum,Csden]=pid_tf(numMs,denMs,mnumT,mdenT,denCs,numCs,denFs);
   %Compute PID from approximated feed-back form controller 
   pid(numMs,denMs,mnumT,mdenT,Csnum,Csden,denCs,numCs);
end % if col > 1

% Disturbance lag------------------
numPd=mt2poly(pdnummt{1,1},x);
denPd=mt2poly(pddenmt{1,1},x);

denCs=conv(denCs,denFs);
Setpoint=1;
Disturb=0;
st_size=1;
Pdelay=pdeadx;
Mdelay=mdeadx;
if ~simulate
   return
end
%-------------------------------------
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_remove(''add'');']);
uicontrol(old_handle,...
  'Style','push',...
  'Position',[555 110 75 25],...
  'String','Change time',...
  'Callback',['add_remove(''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(4260)');
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_1_lag_PID=0; plot_2_lag_PID=0;
 
 IMC_Enable='on'; real_PID_Enable='on'; lag1_PID_Enable='on'; lag2_PID_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,50,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_remove(''add'');']);
uicontrol(f,...
  'Style','push',...
  'Position',[555 110 75 25],...
  'String','Change time',...
  'Callback',['add_remove(''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(4260)');
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);
y_scale=-0.2;
y_position=0.3;
x_position=0.7*time;
options = simset('SrcWorkspace','current');
temp=Epsilon(1);
axis('auto');
if (plot_IMC)
   if plot_setpoint
   	Setpoint=1;
   	Disturb=0;
   	if col > 1
    		if unstable
      	   [t]=sim('imc1dfu2',time,options);
       	  temp1='* IMC (feedback form)'; 
   		else
      		[t]=sim('imc1df2',time,options);
        	 temp1='* IMC'; 
   		end
   	else
      	if unstable
      		[t]=sim('imc1dfu',time,options);
         	temp1='* IMC (feedback form)'; 
   		else
      		[t]=sim('imc1df',time,options);
         	temp1='* IMC'; 
         end
      end
      if plot_effort
      	subplot(2,1,1);plot(t,simout,'m');
        	title(['Setpoint Step Response when Epsilon = ',num2str(temp)])
   		xlabel('Time'); ylabel('Output'); 
         if isempty(keep) | ~keep
            h=text(x_position,y_position,temp1);
            set(h,'color','m');
	   		y_position=y_position-y_scale;
         end
   		grid on;
         hold on
         
        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
      	subplot(2,1,1);plot(t,simout,'m');
        title(['Setpoint Step Response when Epsilon = ',num2str(temp)])
   		xlabel('Time'); ylabel('Output'); 
         if isempty(keep) | ~keep
            h=text(x_position,y_position,temp1);
            set(h,'color','m');
	   		y_position=y_position-y_scale;
         end
   		grid on;
         hold on
      else
         plot(t,simout,'m');
        title(['Setpoint Step Response when Epsilon = ',num2str(temp)])
   		xlabel('Time'); ylabel('Output'); 
         if isempty(keep) | ~keep
            h=text(x_position,y_position,temp1);
            set(h,'color','m');
	   		y_position=y_position-y_scale;
         end
   		grid on;
         hold on
      end
      
   end
   
   if plot_disturb
   	Setpoint=0;
   	Disturb=1;
   	if col > 1
    		if unstable
        	 	[t]=sim('imc1dfu2',time,options);
         	temp1='* IMC (feedback form)'; 
   		else
      		[t]=sim('imc1df2',time,options);
         	temp1='* IMC'; 
   		end
   	else
      	if unstable
      		[t]=sim('imc1dfu',time,options);
         	temp1='* IMC (feedback form)'; 
   		else
      		[t]=sim('imc1df',time,options);
         	temp1='* IMC'; 
      	end
      end
      if plot_effort
      	subplot(2,1,1);plot(t,simout,'m');
        title(['Disturbance Step Response when Epsilon = ',num2str(temp)])
   		xlabel('Time'); ylabel('Output'); 
         if isempty(keep) | ~keep
            h=text(x_position,y_position,temp1);
            set(h,'color','m');
	   		y_position=y_position-y_scale;
         end
   		grid on;
         hold on