# ECE 486 Control Systems Lab (Fall 2017)

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 6^{1} 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}

^{2}