ECE 486 Control Systems Lab (Fall 2017)

[Back to Home]

Day 1 Day 2 Day 3 Day 4 Short labs, weekly
Day 5 Day 6 Day 8   Long labs, biweekly
Day 7 Day 9 Day 10 Day 11 Final project, weekly

Lab 5 – PD Control: Analog Computer and Windows Target

Day 61 of ECE 486 Lab. We continue our experiment with DC motor assembly and further investigate its behavior when we apply different PD controllers.

Maths Behind the Scenes

The name of “PD Control”, means control effort \(u(e(t)) := K_p e(t) + K_d \frac{d }{dt}\{e(t)\}\), where \(K_p\) term is the proportional control and \(K_d\) term is the derivative control. \(u(e(t))\) is a signal depending on error \(e(t) := - y(t) + r(t)\), the difference between actual output and reference input.

Recall from Lab 3, for a second order transfer function (prototype), a high \(K_p\) will suppress disturbance output and by introducing \(K_d\) we can manipulate damping ratio (hence controlling overshoot in reference output).

In this lab, different combinations of \(K_1\) and \(K_2\) (or combinations of \(P_1 = \frac{1}{10} K_1\) and \(P_2 = \frac{1}{10} K_2\)) correspond to different PD controllers, with \(K_1\) being the proportional control gain and \(K_2\) the derivative control gain.

Matlab Part

Keep the header/preamble of your script from Lab 0. (Include those lines as part of your matlab script template.)

In the last section of the lab, you may see the step response will not settle before it changes directions. Therefore for the post lab, you may find the following “projecting” script useful.2 Note it is a function without input arguments. You can use “Publish” button to generate an *.html page as this on the course website.

%% How to estimate $K$, $\zeta$ and $\omega_n$ of a second order TF with partial data?
% estimate the steady state value of the step response of a second order TF
% with partial data

%%

% 2016-10-24
% Y\"un Han
% ECE 486 Lab 5

function parmEstimate2ndOrderTF
%% preamble
clear     % clear values of variables in workspace
clc       % clear messages in the command window
clf       % clear existing figures 
close all % close all existing windows; w/o 'all', only close the latest

%% generate step response with given parameters K, wn and zeta
% K    - DC gain
% wn   - natural frequency
% zeta - damping ratio

% choose some K, wn and zeta
K    =  3.5;
wn   =   10;
zeta = 0.05; % choose low damping so that it takes longer to settle

% we can always generate a unit step response using TF defined by K, wn, zeta
sys = tf(K*wn^2, [1 2*zeta*wn wn^2]);
% obtain the data for step response of the 2nd order TF
Tf = 20; % collect the step response between 0 and Tf
[y, t] = step(sys, 0:0.002:Tf); % UNIT response data
% plot step response
figure(1)
plot(t, y, 'r')
xlabel('time [s]')
ylabel('step response')
title('step response data (t = 0 through 20)')

%% estimate the full step response graph given partial data from figure 1
% the problem in Lab 5 was, as seen in figure 1, the response was highly
% oscillatory and when we collected data, we gave a short time window. say
% in figure 1, we only collected the step response data for a time window 4 sec. 
trimIdx = find(t>4, 1, 'first');
figure(2)
plot(t(1:trimIdx),y(1:trimIdx),'b')
xlabel('time [s]')
ylabel('step response')
title('collected step response data (t = 0 through 4)')
% the question is: is it possible to use the graph between t = 0 and t = 4
% to estimate K, wn and zeta of the 2nd order TF such that we can project
% the step response graph from t = 4 on? if possible, we can concatenate
% the projected graph between t = 4 and t = 20, and the data we already got
% between t = 0 and t = 4; the combination of the two would complete the step
% response of a 2nd order TF, i.e., 
% (partial real data + projection = complete response), then you can call your
% script to calculate Mp, tr, ts.

%% use the data shown in figure 2 as an example, to estimate K, wn and zeta
% sampled five peak values in figure 2 using data cursor
% t1 = 0.314;
% y1 = 6.491;
% t2 = 0.944;
% y2 = 5.683;
% t3 = 1.572;
% y3 = 5.094;
% t4 = 2.202;
% y4 = 4.664;
% t5 = 2.832;
% y5 = 4.35;

t1 = 0.314;
y1 = 6.491;
t2 = 0.944;
y2 = 5.683;
t3 = 1.572;
y3 = 5.094;
t4 = 2.202;
y4 = 4.664;
t5 = 2.832;
y5 = 4.35;
tSample = [t1 t2 t3 t4 t5]';
ySample = [y1 y2 y3 y4 y5]';
% show the sample peaks in figure 2
figure(3)
plot(t(1:trimIdx),y(1:trimIdx),'b')
hold on
plot(tSample,ySample,'linestyle','none', ...
    'marker','s','markersize',8, ...
    'markeredgecolor','red','markerfacecolor',[1 0.5 0.5])
% mark down the peaks
annotation('textarrow',[3.2/4.5   0.2],[6/7 ySample(1)/7.5], 'string',' ')
annotation('textarrow',[3.2/4.5 0.309],[6/7 ySample(1)/8.3], 'string',' ')
annotation('textarrow',[3.2/4.5 0.418],[6/7 ySample(1)/9.1], 'string',' sampled peaks')
annotation('textarrow',[3.2/4.5  0.52],[6/7 ySample(1)/9.7], 'string',' ')
annotation('textarrow',[3.2/4.5  0.63],[6/7 ySample(1)/10.2],'string',' ')
xlabel('time [s]')
ylabel('step response')
title('collected step response data (t = 0 through 4)')

% use (t1, y1) through (t5, y5) to estimate K, wn and zeta
% the time domain solution of unit response of second order TF 
% y(t)=K*(1-1/sqrt(1-zeta^2)*exp(-zeta*wn*t)*sin(wn*sqrt(1-zeta^2)*t+acos(zeta)))
% --- let's say it is eqn.(+)

% denote damped natural frequence as wd = wn*sqrt(1-zeta^2)
% period is the time difference between adjacent peaks, T = 2*pi/wd
Tavg = mean(diff(tSample));
wd = 2*pi/Tavg; % = wn*sqrt(1-zeta^2), so we know wn times some zeta equals a number

% use data points (t1,y1), (t2,y2), (t3,y3)
% denote M = 1/sqrt(1-zeta^2)*exp(-zeta*wn*t1)*sin(wd*t1+acos(zeta))
%        N = exp(-zeta*wn*T)
% by eqn.(+) above, we have
% y1 - y2 = K*M*N - K*M = K*M*(N-1)         --- eqn.(2)
% similarly,
% y1 - y3 = K*M*N^2 - K*M = K*M*(N-1)*(N+1) --- eqn.(3)
% divide eqn.(3) by eqn.(2), we have
% N + 1 = (y1-y3)/(y1-y2)
% so N = (y1-y3)/(y1-y2) - 1                --- eqn.(N cal.)
% note that N by this method, is independent of M, we can repeat this
% using data points (t2,y2), (t3,y3), (t4,y4)
%                   (t3,y3), (t4,y4), (t5,y5) ...
%                   (t(n),y(n)), (t(n+1),y(n+1)), (t(n+2),y(n+2))

N1 = (y1-y3)/(y1-y2) - 1;
N2 = (y2-y4)/(y2-y3) - 1;
N3 = (y3-y5)/(y3-y4) - 1;
Navg = mean([N1 N2 N3]);
% according to the notation of N
% now we have
% zeta*wn = -log(N)/T
% also by notation of wd
% sqrt(1-zeta^2)*wn = 2*pi/T

% solve for wn and zeta by using N and T
% we need to call numerical solver fsolve for a system of nonlinear
% equations
% note: you need Optimisation Toolbox in order to use fsolve
x0 = [.5; 2];  % init guess for zeta and wn 
% call solver; check funcWnZetaVal, it should be almost 0
% parmSolve2ndOrderTF is defined below
if (strcmp(version,'7.14.0.739 (R2012a)'))
    options = optimset('Display','iter'); % for the older matlab 
else
    options = optimoptions('fsolve','Display','iter'); % display iteration output
end
[x,funcWnZetaVal] = fsolve(@(x)parmSolve2ndOrderTF(x,Navg,Tavg),x0,options) 


zetaEst = x(1) % compare our estimation with the acutal zeta = 0.05
wnEst   = x(2) % compare with the actual wn = 10

% how about K, the steady state value - the most important parameter we
% care about? K = y1/(1-M) etc
K1 = y1/(1-MCal(zetaEst,wnEst,t1)); % Mcal is defined below
K2 = y2/(1-MCal(zetaEst,wnEst,t2));
K3 = y3/(1-MCal(zetaEst,wnEst,t3));
KEst = mean([K1 K2 K3]) % compare with the actual K = 3.5

% estimation success!

%% complete the partial plot with our projection
figure(4)
% create the projected TF for from t = 4 sec onwards
sysProj = tf(KEst*wnEst^2, [1 2*zetaEst*wnEst wnEst^2]);
% projected response data
[yProj, tProj] = step(sysProj, 0:0.002:Tf); % UNIT response data
% plot both collected data (blue) and projected data (red)
plot(t(1:trimIdx),y(1:trimIdx),'b','linewidth',1.5)
hold on
plot(tProj(trimIdx:end),yProj(trimIdx:end),'r--','linewidth',2)
xlabel('time [s]')
ylabel('step response')
legend('collected data','projected data','location','best')
text(10,2,'is this figure quite similar to figure 1?')
title('collected step response data (t = 0 through 4) and projected data (t = 4 onwards)')
end % parmEstimate2ndOrderTF

function funcWnZeta = parmSolve2ndOrderTF(x,N,T)
    % x = (zeta, wn)
    funcWnZeta = [x(1)*x(2) + log(N)/T;
                  sqrt(1-x(1)^2)*x(2) - 2*pi/T];
end

function M = MCal(zeta,wn,t)
    M = 1/sqrt(1-zeta^2)*exp(-zeta*wn*t)*sin(wn*sqrt(1-zeta^2)*t+acos(zeta));
end

Visualize data using plot().

Hardware Part

  • Since we are going to use Simulink™ Windows Target®, it is not a good idea to run your files from a network path. For example, you may need to create a folder named after your NetID in local C:\ drive and save your Matlab WinTarget files there. Move the folder to your thumb drive once you are finished.

    Always keep C:\ drive clean and organized.

  • When you collect data for step responses of the closed loop system, you can activate “Quick Measure” on the oscilloscope panel and record \(M_p\) and \(t_r\) (no \(t_s\) however!) at the same time. This is a “cheating” way to complete several data tables in the post lab. There, you are supposed to call StepResponseMetrics.m you wrote in Lab 0 to calculate those specs. But you can also use the specs you record from oscilloscope to sanity check.
  • Always clean up before exiting. Specifically,
    • Clean up bench table, restore pot, motor lock etc, reinstall screws;
    • Sort out wires color by color, type by type and put them back to racks;
    • Turn off oscilloscope, meters etc;
    • Restore chairs.

Follow-up

Comments

  • For our next lab (during Week of Oct 30), no prelab or lab report will be due in lab, but all (both Group A and B) students attend lab for Week 1 Final Project. We will do Lab 6 after Week 1 Final Project. The point of doing this is we want the lecture to stay a little bit ahead of us when it comes to frequency domain design. And good news is that prelab for Final Project is due in Week 2 Final Project. Please check course website for the exact date, it should be around Nov 27.
  • When you fill out the tables in post lab, if your step responses took more than half of the period to settle, you may ask yourself the question “how can I get the steady state value?” It is OK to leave a note in your lab report saying that you didn’t see the steady state value directly from the data set. But we actually can estimate the steady state value based on an “incomplete” data set. The technique is illustrated above in the Matlab part. It is optional.
  • Some of you might have seen “wild” motor responses when full compensation was applied. Note that the wild response is stable, since position reading is bounded by \(\pm 10\) Volt. The “wildness” was caused by modular arithmetics. The potentiometer voltage reading is modded by 10, with \(0\) Volt corresponding to \(0^\circ\), \(10^-\) Volt corresponding to \(360^{-\circ}\). After friction compensation is applied, the motor is so responsive in some cases that large overshoot is incurred initially. If the peak position voltage exceeds \(10\) Volt, say we have a positive rotation of \(540^\circ\), the potentiometer shall give \(15\) Volt, but it only returns \(15 \bmod 10 = 5\) Volt. Therefore controller will go ahead and carry out calculation of the desired control effort based on this reading \(5\) Volt but not \(15\) Volt, i.e., the controller is blind to modular arithmetics. From this point on, the control calculation has lost track. This causes the wild behavior.

Due Date

Lab 5 Report is due at the beginning of Lab 6 (Nov 9, Group A; Nov 16, Group B). No prelab is due on Monday by 5pm of the week of Week 1 Final Project. Prelab 6 is due on Monday (Nov 6, Group A; Nov 13, Group B) in the week when you have Lab 6. Handwritten prelab is acceptable; typesetting is recommended.

Questions

You are always very welcome to stop by office hours on Mondays. Emailing questions is another way. You can always include [ECE486]blah in the title of your question emails.

Spot any typos? Email me at once. You will earn up to +5 points for each typo/technical error reported.

Footnotes:

1
Lab class meets every other week until final project starts.
2
It is useful in the sense that you may use it for completion of Table 4 of Lab 5 Report.

Author: Yün Han

Last updated: 2017-10-26 Thu 21:12