Course materials · Control Engineering I

Control Engineering I · Week 14 · Comparing Stability and Performance Indices

English translation

Edit the code and response fields in your own copy. Name: ⟦enter here⟧ / Student ID: ⟦enter here⟧

How to run

In Live Editor, click inside the code and run the current section with Run Section (Windows/Linux Ctrl+Enter, Mac ⌘+Enter). Run executes the entire document. Edit the explanations and response fields directly.
Initially run from the top, and rerun any section whose values you change. If an undefined-variable error occurs, run the preceding sections first. Outputs appear below or to the right of the code.
Environment: MATLAB Online and Control System Toolbox.

1. Calculate the closed-loop model and final value

Perform the calculations below first, then compare with MATLAB results.
  • Use G(s)=1s+2G(s)=\frac1{s+2} and H(s)=2s+4H(s)=\frac2{s+4} from the “3. Tracking error and error at the comparison point” code in the practice Live Script, but first calculate the case K=2K=2 by hand.
  • Substitute U=K(R−HY)U=K(R-HY) into Y=GUY=GU. Collect the Y terms on one side to simplify YR\frac{Y}{R}, then find the denominator roots and determine whether all poles have negative real parts.
  • Find the final output and sensor measurement for a unit step, then calculate the final values of r−yr-y and r−Hyr-Hy separately.
  • Then set K to 2, run the code, and compare the two results.
  • Skip the already calculated K=2K=2 and continue with K=4,8K=4,8.
  • For the Routh activity, also derive the first column yourself as previously instructed, then check using roots.
Before execution — Hand-calculated expressions and values: ⟦write calculations here⟧ After execution — MATLAB output values to compare: ⟦write output values here⟧

2. Determine stability

A pole with real part 0 does not meet the condition that every pole has a negative real part. Check where the poles lie at the boundary values.
Prediction before running: ⟦enter here⟧
polynomials={[1 3 2],[1 1 2 8],[1 0 1]};
for j=1:numel(polynomials)
p=roots(polynomials{j});strictStable=all(real(p)<-1e-9);
disp(p);disp(strictStable);
end
Do it yourself:
  1. Derive the first Routh column of s3+3s2+2s+Ks^3+3s^2+2s+K and compare K∈{1,6,7}K\in\{1,6,7\}.
  2. Find the poles of K=6K=6 and check whether any have real part 0.
  3. Explain how the response differs from the case where all poles have negative real parts.
% ✍ Your code: add your own code below.
Changed condition: ⟦enter here⟧
Observed result: ⟦enter here⟧
Differences from the prediction and their causes: ⟦enter here⟧

3. Tracking error and error at the comparison point

y is the actual output, and Hy is the sensor measurement. Subtracting each from the reference gives r−y and r−Hy. Calculate both side by side to check whether the sensor causes a difference.
Prediction before running: ⟦enter here⟧
s=tf('s');K=4;G=1/(s+2);H=2/(s+4);T=feedback(K*G,H);
assert(isstable(T));
tracking=1-dcgain(T);comparator=1-dcgain(H*T);disp([tracking comparator]);
Do it yourself:
  1. For K∈{2,4,8}K\in\{2,4,8\}, calculate r−yr-y, the target minus the actual output, and r−Hyr-Hy, the target minus the sensor measurement, separately.
  2. State which error to use to judge whether the actual output is close to the reference.
% ✍ Your code: add your own code below.
Changed condition: ⟦enter here⟧
Observed result: ⟦enter here⟧
Differences from the prediction and their causes: ⟦enter here⟧

4. Finite-horizon performance index

Integrate with trapz using the same time range and spacing.
Prediction before running: ⟦enter here⟧
s=tf('s');G=1/(s^2+2*s+1);t=(0:.01:20)';
for K=[1 2 4]
T=feedback(K*G,1);[y,~]=step(T,t);e=1-y;
disp([K trapz(t,e.^2) trapz(t,abs(e)) trapz(t,t.*abs(e))]);
end
Do it yourself:
  1. Using the plant G(s)=1/(s2+2s+1)G(s)=1/(s^2+2s+1) and K from the ISE calculation, calculate ISE with end times of 10 s, 20 s, and 40 s.
  2. Tabulate the K used, integration end time, and ISE. For each K, compare how much ISE increases from 10→20 s and from 20→40 s.
  3. Check whether the increments decrease and judge whether the integral approaches a constant. Include the integration interval in your results, and distinguish finite-horizon integrals from infinite-horizon indices.
% ✍ Your code: add your own code below.
Changed condition: ⟦enter here⟧
Observed result: ⟦enter here⟧
Differences from the prediction and their causes: ⟦enter here⟧

5. Optimization and Approximation

State the design range and comparison interval.
Prediction before running: ⟦enter here⟧
K3=.1:.01:10;ISE=1./K3+.1*K3;[v,idx]=min(ISE);disp([K3(idx),sqrt(10),v]);
s=tf('s');GH=1/(1+11*s/6+s^2+s^3/6);GL=1/(1+1.615*s+.624*s^2);
t=(0:.01:15)';yh=step(GH,t);yl=step(GL,t);
figure;plot(t,yh,t,yl);legend('GH','GL');grid on;disp(max(abs(yh-yl)));
Do it yourself:
  1. Keep the K3 search range at 0.1–10, change only the spacing, and rerun the minimum search. Tabulate the search spacing before and after the change, the gain giving the minimum, and the absolute difference between 10\sqrt{10} and that gain.
  2. Create a comparison model by multiplying only one denominator coefficient of GL by 1.1. For example, extract the coefficients using [numL,denL]=tfdata(GL,'v'), apply denL(end)=1.1*denL(end), and create a new model using tf(numL,denL). State the power associated with the changed coefficient and its values before and after the change.
  3. Calculate max(abs(yh-yl)) before and after the change at the same t=0:0.01:15. Use this maximum absolute difference to state whether the approximation error for GH increases or decreases.
% ✍ Your code: add your own code below.
Changed condition: ⟦enter here⟧
Observed result: ⟦enter here⟧
Differences from the prediction and their causes: ⟦enter here⟧

Practice through repetition and check before submission

Summarize the three conditions compared in the preceding hands-on activities below. You do not need to repeat experiments already completed. Predictions may state increasing or decreasing trends instead of exact values.
Repetition 1 — Condition: ⟦enter here⟧ / Prediction: ⟦enter here⟧ / Actual result: ⟦enter here⟧
Trial 2 — Condition: ⟦enter here⟧ / Prediction: ⟦enter here⟧ / Observation: ⟦enter here⟧
Trial 3 — Condition: ⟦enter here⟧ / Prediction: ⟦enter here⟧ / Observation: ⟦enter here⟧
Optional:
Example: Classifying an imaginary-axis pole as asymptotically stable
  1. Find the poles using p=roots([1 0 1]), and record the result of applying all(real(p)<=0) as the decision criterion.
  2. Write this condition beside the asymptotic stability test all(real(p)<-1e-9). Distinguish poles with real part 0 from those with negative real parts, and restore only the conditional expression to test asymptotic stability.
Include code, actual outputs and plots, and your own interpretation. Submit through KLAS; check the semester schedule and latest announcements for the deadline.
Filename: Name_StudentID_Week14.mlx.