function step2demo() % STEP2DEMO -- Demonstrate relationships between pole locations and step response for a second % order system. % Click the mouse on the S plane to locate a complex pole. Its conjugate appears automatically, % and the step response is drawn. Click outside the graph to quit. % % If the pole is close enough to the real axis, it is assumed real, and you are asked to % click a second pole. Critical damping is supported, but not poles at the origin. % % Pole locations are coarsely quantized, for simplicity. % E.W. Hansen, Engs 22, 97W % Modified C.R. Sullivan 7/22/98 -- direct calculation of response. % Modified E.W. Hansen 1/30/99 -- % Added performance measures, lines of constant wn and Mp in s plane % Modified 7/99 CR Sullivan Added manual/auto scale option, and intro dialog box. button1=strcat('Fixed', '(for easy comparisons)') button2=strcat('Variable', '(for best detailed view of a particular plot)') ascale = menu('Choose a plot scaling option:',button1,button2); ascale = ascale-1 % 0 for fixed, 1 for auto menu(strvcat('An s-plane plot will be displayed.',... 'You can then click on the s-plane to place poles.',... 'To quit, click off the graph (but in the figure window).'),'OK'); % Initialize MaxSigma = 3; MaxOmega = 2; PlotWidth = 2*MaxOmega; corner = sqrt(MaxSigma^2 + MaxOmega^2); Nt = 200; th = linspace(pi/2, 3*pi/2, 51); % Lines of constant wn semicirc = exp(j*th); Mp = [0.01, 0.1, 0.25]; % Lines of constant Mp zeta = abs(log(Mp)) ./ sqrt(pi^2 + log(Mp).^2); q = 0.05; % Quantization on pole locations OD = 1; % Overdamped, critically damped, underdamped CD = 2; UD = 3; h0 = []; % Initialize object handles h1 = []; h2 = []; A = [-MaxSigma, -MaxSigma+PlotWidth, -MaxOmega, MaxOmega]; clf subplot(1,2,1), plot([-MaxSigma, MaxSigma/2], [0, 0], 'm', [0, 0], [-MaxOmega, MaxOmega], 'm'); grid hold on for wn = [1, 2, 3], plot(wn*semicirc, 'k--'); % Constant wn lines end for z = zeta, th = pi-acos(z); plot([0, corner*cos(th)], [0, corner*sin(th)], 'k--') end hold off axis(A), axis('square') title('CLICK ON A POLE LOCATION') % Big event loop while(1) % GET POLE LOCATIONS OFF THE GRAPH subplot(1,2,1) [s, w] = ginput(1); % Get a pole location off the graph s = q * round(s/q); % Round the locations w = q * round(w/q); if (sA(2) | wA(4)) % Outside the frame? name = inputdlg('Enter text to add to the display: '); if ~isempty(name), text(-2.7,-1.8,name); end return end if exist('h1'), set(h1, 'Visible', 'off'); end % Clear previous marks if exist('h2'), set(h2, 'Visible', 'off'); end if w==0, % Real pole? hold on h0 = plot(s, w, 'kx'); hold off; title('SELECT A SECOND REAL POLE') [s2, w2] = ginput(1); s2 = q * round(s2/q); if s==s2, damping = CD; else damping = OD; end s = [s, s2]; w = [0, 0]; polelocation = [num2str(s(1)), ' , ', num2str(s(2))]; else damping = UD; s = [s, s]; % Make a conjugate pair w = [w, -w]; polelocation = [num2str(s(1)), ' ± j',num2str(abs(w(1)))]; end p=s+j*w; % Vector of pole locations a=-p(1); b=-p(2); wn=sqrt(a*b); % Undamped natural frequency zeta=(a+b)/2/wn; % Damping ratio if zeta<1, polelocation = [polelocation, '; \zeta = ', num2str(zeta,3), ', \omega_n = ', num2str(wn,3)]; end if exist('h0'), set(h0, 'Visible', 'off'); end hold on h1 = plot(s, w, 'kx'); title(['s = ', polelocation]); hold off if(damping==UD), hold on, h2 = plot([s(1), 0, s(2)], [w(1), 0, w(2)], 'r'); hold off end xlabel('Real(s)'), ylabel('Imag(s)') % COMPUTE AND PLOT THE RESPONSE if(min(abs(s))==0 & min(abs(w))==0), % System poles at the origin? title('SELECT SYSTEM POLES NOT AT THE ORIGIN') else % p=s+j*w; % a=-p(1); % b=-p(2); % wn=sqrt(a*b); % zeta=(a+b)/2/wn; if ascale if zeta~=0, % Poles off imaginary axis tmax = 8/min(abs(s)); % 8 time constants for settling else tmax = 8*2*pi/wn; % or 8 full cycles end else tmax = 15; %fixed scale end t = linspace(0, round(tmax), Nt); % Time vector for plotting if zeta~=1, % If we're not critically damped % Response from the table #10. ystep =real((1 - b/(b-a)*exp(-a*t) + a /(b-a)*exp(-b*t))) ; % Response from table #2 yimp =real( a*b/(b-a) * (exp(-a*t) - exp(-b*t)) ); else ystep = 1 - (1 + wn*t).*exp(-wn*t); % Critically damped response yimp = wn^2*t.* exp(-wn*t); end % Compute performance measures, where appropriate if zeta<1 [ypeak, ipeak] = max(ystep); Mp = ypeak-1; % Peak overshoot tp = t(ipeak); % Peak time wd = abs(w(1)); % Ringing frequency else Mp = 0; end n10 = max(find(ystep(1:ipeak)<=0.1)); % Risetime n90 = max(find(ystep(1:ipeak)<= 0.9)); tr = t(n90) - t(n10); ns = max(find(abs(ystep-1) >= 0.01)); % Settling time ts = t(ns); subplot(1,2,2), plot(t, ystep, 'b', t, yimp, 'g'), grid axis('square') if not(ascale), axis([0,15, -0.2, 1.8]), end xlabel('time'), ylabel('responses') if(damping==OD), titlestring = ['OVERDAMPED, t_r=', num2str(tr, 3)]; else if (damping == CD), titlestring = ['CRITICALLY DAMPED, t_r=', num2str(tr, 3)]; else % titlestring = ['\zeta = ', num2str(zeta,3), ', \omega_n = ', num2str(wn,3)]; if Mp>0, titlestring = ['PO=',num2str(Mp,2), ', ']; else titlestring = []; end titlestring = [titlestring, 't_r=',num2str(tr,3), ', t_s=',num2str(ts,3),... ', \omega_d=', num2str(wd,3)]; end end title(titlestring) end end