ex8_3a.m

离散控制系统设计的MATLAB 代码

%%%%%%%%%%%%%%%%% Example 8.3(a) %%%%%%%%%%%%%%%%
%   Discrete-Time Control Problems using        %
%       MATLAB and the Control System Toolbox   %
%   by J.H. Chow, D.K. Frederick, & N.W. Chbat  %
%         Brooks/Cole Publishing Company        %
%                September 2002                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   ----  Proportional-integral-derivative control  ----
%
clear
disp('Example 8.3(a)')

Ts = 0.1;			                % sampling period
Gp_s = tf(4,conv([2 1],[0.5 1]));   % continuous process
H_s = tf(1,[0.05 1]);		        % continuous sensor
Gp_z = c2d(Gp_s,Ts,'zoh');	        % discrete plant (no sensor)
HGp_s = H_s*Gp_s                    % continous plant
HGp_z = c2d(HGp_s,Ts,'zoh')         % discrete plant, using ZOH option

dtime  = 0:Ts:2;
KI = 0.5; 	
symbols = ['o' '+' '*' 'd' 'v'];    % different plot symbols
for KD = 0.3:0.1:0.6
  figure
  results = [];
  Pgains = [];
  hold on
  for ii = 1:4 		                % KP values are: 1.7, 1.9, 2.1, 2.3
    KP = 1.5+0.2*ii;
    %-- write Gc(z) as sum of 3 LTI objects in TF form --
    Gc = KP*(tf(1,1,Ts) ...			% porportional
	+ tf([KI*Ts  0],[1 -1],Ts) ...	% integral
	+ tf(KD*[1 -1],[Ts 0],Ts)); 	% derivative

    % reference input - CL transfer fcn
    Tref_z = Gp_z*Gc/(1+HGp_z*Gc);
    Tref_z = minreal(Tref_z);     % remove pole-zero cancellations
    y = step(Tref_z,dtime);       % CL response to step input
    ref = dcgain(Tref_z);
    [Mo,tpeak,trise,tsettle,ess] = kstats(dtime,y,ref);
    ymax = y(find(dtime == tpeak));
    yss = y(length(dtime));
    results = [results; [KP KD trise tsettle  Mo]];
    Pgains = [Pgains KP];
    disp(['KP = ',num2str(KP)])
    plot(dtime,y,symbols(ii))
    pause(1)
  end  % for ii..
  hold off                          % turn hold off after last plot
  xlabel('Time (s)')
  ylabel('Step Response')
  title(['Ex8.3 a - Step Responses for Different ',...
        'Values for KP, where KI = ',num2str(KI),...
        ' and KD = ',num2str(KD)]) 
  grid on
  text('Units','norm','Pos',[0.2,0.15],'Str',['KI = ',num2str(KI)])
  text('Units','norm','Pos',[0.2,0.10],'Str',['KD = ',num2str(KD)])

  legend('KP = 1.7','KP = 1.9','KP = 2.1','KP = 2.3')

  %----- print results table
  disp(['KI = ',num2str(KI),' and KD = ',num2str(KD)]) 
  disp(' ')
  disp('      KP       KD     t_rise      tsettle     M_o (pct)')
  disp('    ---------------------------------------------------')
  disp(results)
  disp('    ---------------------------------------------------')
  disp(' ')
end  % for KD..
%%%%%%%%%%