Root Locus Plot Is Way Off

27 views (last 30 days)
Arpad
Arpad on 26 Nov 2023
Edited: Paul on 2 Dec 2023
If anyone knows how to use matlab, can someone explain what is going wrong.
I am working on a project for a class called Linear Systems. I am plotting a root locus and it is supposed to have a curve that goes through 2 desired poles but its so off and its not even funny. I have a transfer function with desired poles.
This is my code: (The first section is very similar to my professors. Its the second section that is just way off).
sys = tf(0.03965,[5.28e-08,1.856e-05,0.00080187,0]);
rlocus(sys)
hold on;
plot(-7.9,11.85,'rx')
plot(-7.9,-11.85,'rx')
% sys_new = tf([10],[1 5])*tf([1],[1 60])*tf([1],[1 0])*tf([1 13.33],[1])
sys_new = tf(0.03965,[5.28e-08, 1.856e-05, 0.00080187, 0])*tf([1 22.835],1);
figure()
rlocus(sys_new)
hold on;
% plot(-20,20,'rx')
% plot(-20,-20,'rx')
plot(-7.9,11.85,'rx')
plot(-7.9,-11.85,'rx')
The part of the code with % is my professors.

Answers (4)

Paul
Paul on 26 Nov 2023
Hi Arpad,
Let's take a look at that root locus for the uncompensated system. I asusme that the red x's indicate the desired closed loop pole locations after compensation.
sys = tf(0.03965,[5.28e-08,1.856e-05,0.00080187,0]);
rlocus(sys)
hold on;
plot(-7.9,11.85,'rx')
plot(-7.9,-11.85,'rx')
Now add the compensation, which is a single zero at s = -22.835. Also, it's better use sys to define sysnew to avoid typos or copy/paste errors.
% sys_new = tf([10],[1 5])*tf([1],[1 60])*tf([1],[1 0])*tf([1 13.33],[1])
%sys_new = tf(0.03965,[5.28e-08, 1.856e-05, 0.00080187, 0])*tf([1 22.835],1);
sys_new = sys*tf([1 22.835],1);
figure()
rlocus(sys_new)
hold on;
% plot(-20,20,'rx')
% plot(-20,-20,'rx')
plot(-7.9,11.85,'rx')
plot(-7.9,-11.85,'rx')
So I guess the question is what is the basis for thinking that adding the zero at s = -22 would cause the root locus to go through, or even near, the desired pole locations. As shown, the compensator zero results in a real closed loop pole somewhere between the zero and the plant pole near the orign. If you put the compensator zero to the left of the plant pole at s = -50, the root locus will qualitatively look more like the "proferssor's solution", but the loci will be pulled much further to the left than the desired closed-loop poles. So maybe a single compensator zero isn't sufficient for this plant and the desired closed-loop pole locations.
  26 Comments
Paul
Paul on 27 Nov 2023
In your original code you had this term
tf([1 22.835],1)
ans = s + 22.84 Continuous-time transfer function.
which has one zero at s0 = -22.835 and no poles. So you need to figure out how to change that tf command to represent a pole at s0 = -22.835 and no zeros. If you don't know how to that, I suggest checking the doc page for tf. Or, it might be more intuitive to use a zpk form and convert it to a tf. For example
tf(zpk(-22.835,[],1)) % same as above.
ans = s + 22.84 Continuous-time transfer function.
Strictly speaking, it's not necssary to convert the zpk to a tf. I just did that here for comparison to your original tf.
BTW, that was a nice idea to use the angle criterion to determine the pole location of the compensator.
Paul
Paul on 2 Dec 2023
Edited: Paul on 2 Dec 2023
The Symbolic Toolbox can be helpful in solving a problem like this.
Define the plant and the desired location of the dominant pair of complex poles.
sys = tf(0.03965,[5.28e-08,1.856e-05,0.00080187,0]);
pdes = -7.9+1j*11.85;
Define the plant symbolically
syms s
P(s) = simplify(poly2sym(sys.num{1},s)/poly2sym(sys.den{1},s));
Asusme we want to us a single-pole compensator
syms p real
K(s) = 1/(s + p);
Solve for the value of p that satifsfies the root locus angle criteria
sol.p = solve(angle(P(pdes)*K(pdes))==sym(pi),p);
double(sol.p)
ans = 22.8346
Same value as found above. Plot the root locus using this compensation
figure
rlocus(sys*tf(1,[1 double(sol.p)])),
hold on
plot(pdes,'r*')
plot(conj(pdes),'r*')
axis([-30 0 -70 70])
The compensator pole introduces a breakaway point on the real axis to the left of the desired poles and the loci migrate to the right through the desired pole locations as the gain increases.
Assume we want a compensator with one zero and one pole
syms z real
K(s) = (s+z)/(s+p);
Solve for the locations based on the angle criteria. Because we have one equation in two unknowns, we solve for one in terms of the other. I chose to solve for p in terms of z.
sol = solve(angle(P(pdes)*K(pdes))==sym(pi),p,'ReturnConditions',true);
Pick a value of z. Here, I'm using the value of z as was originally defined in the question.
z = 22.835;
Get the value of p that corresponds to that value of p.
p = double(subs(sol.p,z))
p = 10.6662
Plot the root locus.
figure
rlocus(sys*tf([1 z],[1 p])),
hold on
plot(pdes,'r*')
plot(conj(pdes),'r*')
axis([-30 0 -70 70])
We see that shifting the compensator pole to lower frequency yields a breakaway point the right of the desired pole locations and the compensator plays its role as an attractor to bring the loci to the left through the desired pole locations before the loci evetually migrate toward the RHP as the gain increases.

Sign in to comment.


Sam Chak
Sam Chak on 26 Nov 2023
I am uncertain about what your professor intends for you to learn. Initially, the problem seems akin to a standard Pole Placement issue. However, upon delving into it for several hours, it emerges as a mathematically challenging problem, particularly if you aim to precisely align the root locus to traverse the two specified closed-loop poles while exhibiting the desired 2nd-order system behavior simultaneously.
From the given complex poles , a reference 2nd-order system can be formulated as . Although you suggested an improper compensator, for practical systems, a strictly proper compensator is preferable, wherein the numerator's degree is lower than that of the denominator.
Essentially, as previously mentioned, this problem can be tackled using pole placement only when all desired poles are known, as is the case in most exercises. However, the plant system is of 3rd-order, and only two desired closed-loop poles are known. Hence, the dominant pole approximation method can be employed. You can explore this topic on Google, or it may already be covered in your professor's lectures.
In my proposed solution, I opted for a relatively high-order but strictly proper compensator.
s = tf('s');
% 3rd-order Plant transfer function
Gp = tf(0.03965, [5.28e-8 1.856e-5 0.00080187 0])
Gp = 0.03965 ------------------------------------------ 5.28e-08 s^3 + 1.856e-05 s^2 + 0.0008019 s Continuous-time transfer function.
% A relatively high-order Compensator (but strictly proper)
a1 = 31264.2848484848;
a2 = 363954426.075654;
a3 = 1849668573733.87;
b1 = 4355207983.5442;
b2 = 45076892441.6564;
b3 = 1052053220570.43;
c1 = 241.562107836302;
d1 = 10.3501124658055;
d2 = c1;
G1 = (b1*s^2 + b2*s + b3)/(s^3 + a1*s^2 + a2*s + a3);
G2 = tf(c1, [1 d1 d2]);
% My compensator formula (maybe it appears in literature, I didn't check)
Gcom = G2*G1/(1 + G1*Gp - G2*G1*Gp);
Gcom = minreal(Gcom)
Gcom = 1.052e12 s^14 + 6.654e16 s^13 + 1.842e21 s^12 + 2.914e25 s^11 + 2.815e29 s^10 + 1.61e33 s^9 + 4.669e36 s^8 + 2.848e39 s^7 + 6.29e41 s^6 + 5.219e43 s^5 + 1.982e45 s^4 + 4.307e46 s^3 + 7.17e47 s^2 + 6.413e48 s + 4.844e49 ----------------------------------------------------------------------------------------------------------------------------------------------------------------------- s^17 + 9.452e04 s^16 + 4.092e09 s^15 + 1.073e14 s^14 + 1.891e18 s^13 + 2.342e22 s^12 + 2.066e26 s^11 + 1.282e30 s^10 + 5.376e33 s^9 + 1.393e37 s^8 + 1.812e40 s^7 + 5.524e42 s^6 + 3.969e44 s^5 + 1.214e46 s^4 + 2.585e47 s^3 + 3.599e48 s^2 + 3.242e49 s + 1.878e50 Continuous-time transfer function.
% Closed-loop system
Gcls = minreal(feedback(Gcom*Gp, 1))
Gcls = 7.9e17 s^14 + 4.997e22 s^13 + 1.383e27 s^12 + 2.189e31 s^11 + 2.114e35 s^10 + 1.209e39 s^9 + 3.506e42 s^8 + 2.139e45 s^7 + 4.723e47 s^6 + 3.919e49 s^5 + 1.489e51 s^4 + 3.234e52 s^3 + 5.384e53 s^2 + 4.816e54 s + 3.638e55 ----------------------------------------------------------------------------------------------------------------------------------------------------------------------- s^20 + 9.487e04 s^19 + 4.126e09 s^18 + 1.087e14 s^17 + 1.929e18 s^16 + 2.409e22 s^15 + 2.148e26 s^14 + 1.355e30 s^13 + 5.83e33 s^12 + 1.584e37 s^11 + 2.31e40 s^10 + 1.211e43 s^9 + 2.617e45 s^8 + 2.377e47 s^7 + 1.102e49 s^6 + 3.18e50 s^5 + 6.712e51 s^4 + 9.858e52 s^3 + 1.097e54 s^2 + 7.668e54 s + 3.638e55 Continuous-time transfer function.
% Check closed-loop system poles
CL_poles = pole(Gcls);
CL_poles(19:20)
ans =
-7.9000 +11.8500i -7.9000 -11.8500i
% Plot Step response
figure(1)
step(Gcls, 1), hold on
% Compare with the reference 2nd-order Transfer function
Gref = 202.833/(s^2 + 15.8*s + 202.833)
Gref = 202.8 -------------------- s^2 + 15.8 s + 202.8 Continuous-time transfer function.
step(Gref, 1), grid on, hold off
legend('Closed-loop system', '2nd-order system')
% Plot Root Locus
figure(2)
rlocus(Gcom*Gp), hold on
plot(-7.9, 11.85, 'rx', 'markersize', 12, 'linewidth', 3)
plot(-7.9, -11.85, 'rx', 'markersize', 12, 'linewidth', 3)
axis([-10 -6 -20 20]), hold off

Sam Chak
Sam Chak on 26 Nov 2023
It seems that employing the place() command is a more efficient approach to solving the pole placement problem. I will leave the task of converting the matrix gain K to the compensator transfer function as an exercise for you. Please verify this with your professor.
% Plant
Gp = tf(0.03965, [5.28e-8 1.856e-5 0.00080187 0])
Gp = 0.03965 ------------------------------------------ 5.28e-08 s^3 + 1.856e-05 s^2 + 0.0008019 s Continuous-time transfer function.
% Convert Plant TF to State-Space
sys = ss(Gp)
sys = A = x1 x2 x3 x1 -351.5 -118.6 0 x2 128 0 0 x3 0 1 0 B = u1 x1 64 x2 0 x3 0 C = x1 x2 x3 y1 0 0 91.67 D = u1 y1 0 Continuous-time state-space model.
% Make desired closed-loop poles
coef = conv([1 15.8 202.833], [1 7.9e3]); % dominant pole technique
p = roots(coef)
p =
1.0e+03 * -7.9000 + 0.0000i -0.0079 + 0.0119i -0.0079 - 0.0119i
% Find gain matrix K to place closed-loop poles at desired locations
A = sys.A;
B = sys.B;
K = place(A, B, p)
K = 1×3
118.1920 13.4077 195.6031
% Closed-loop system
Acl = A - B*K;
syscl = ss(Acl, B, sys.C, sys.D);
N = 1/dcgain(syscl); % Input scaler
syscl = ss(Acl, B*N, sys.C, sys.D);
Pcl = pole(syscl) % if K is correctly designed, then Pcl = p
Pcl =
1.0e+03 * -7.9000 + 0.0000i -0.0079 + 0.0119i -0.0079 - 0.0119i
% Plot Step response
step(syscl), hold on
% Compare with the reference 2nd-order Transfer function
s = tf('s');
Gref = 202.833/(s^2 + 15.8*s + 202.833)
Gref = 202.8 -------------------- s^2 + 15.8 s + 202.8 Continuous-time transfer function.
step(Gref, 1), grid on, hold off
legend('Closed-loop system', '2nd-order system')

Sam Chak
Sam Chak on 26 Nov 2023
Based on the description from the OP, the compensated system should form a closed-loop system. The formula for a closed-loop transfer function, when the compensator is placed in series with the plant (ys), is given by:
% closed-loop transfer function (CLTF):
% compensated_sys = (compensator*ys)/(1 + compensator*ys) % actual G/(1 + G)
% desired_sys = (compensator*ys)/(1 + compensator*ys) % desired G/(1 + G)
The OP aims to design the compensator in a way that two of the closed-loop poles are (-7.9 ± 11.85i). When these poles are plotted on the root locus diagram, the locus should traverse the two desired closed-loop poles.
s = tf('s');
ys = tf(0.03965, [5.28e-08, 1.856e-05, 0.00080187, 0])
ys = 0.03965 ------------------------------------------ 5.28e-08 s^3 + 1.856e-05 s^2 + 0.0008019 s Continuous-time transfer function.
% Desired system (closed-loop) to maintain unity (1) when s = 1
desired_sys = ((-7.9 - 11.85i)*(-7.9 + 11.85i))/((s - (-7.9 - 11.85i)) * (s - (-7.9 + 11.85i)))
desired_sys = 202.8 -------------------- s^2 + 15.8 s + 202.8 Continuous-time transfer function.
% Compensator (improper, numerator's degree > denominator's degree
% The compensator can be inversely obtained by solving the CLTF equation
compensator = desired_sys/(ys - ys*desired_sys)
compensator = 5.655e-13 s^8 + 4.065e-10 s^7 + 9.344e-08 s^6 + 7.493e-06 s^5 + 0.0002435 s^4 + 0.003285 s^3 + 0.02645 s^2 ---------------------------------------------------------------------------------------------------------- 2.094e-09 s^7 + 8.021e-07 s^6 + 5.6e-05 s^5 + 0.001344 s^4 + 0.01674 s^3 + 0.1019 s^2 Continuous-time transfer function.
% Closed-loop system
compensated_sys = feedback(compensator*ys, 1); % can use feedback() command for CLTF
compensated_sys = minreal(compensated_sys)
compensated_sys = 202.8 s^6 + 1.458e05 s^5 + 3.352e07 s^4 + 2.688e09 s^3 + 8.733e10 s^2 + 1.178e12 s + 9.489e12 ----------------------------------------------------------------------------------------------------------------- s^8 + 734.6 s^7 + 1.768e05 s^6 + 1.601e07 s^5 + 6.735e08 s^4 + 1.53e10 s^3 + 2.259e11 s^2 + 1.918e12 s + 9.489e12 Continuous-time transfer function.
The poles (or roots of the denominator) of the compensated system can be checked using this command:
% Check closed-loop poles
CL_poles = pole(compensated_sys)
CL_poles =
1.0e+02 * -3.0107 + 0.0000i -3.0107 + 0.0000i -0.5044 + 0.0000i -0.5044 + 0.0000i -0.0790 + 0.1185i -0.0790 - 0.1185i -0.0790 + 0.1185i -0.0790 - 0.1185i
% Plot the step response of Compensated system and compare with Desired system
figure(1)
step(compensated_sys, 1), hold on
step(desired_sys, 1), grid on, hold off
legend('Compensated system', 'Desired system')
% Plot Root locus
figure(2)
rlocus(compensator*ys), hold on
plot(-7.9, 11.85, 'rx', 'markersize', 12, 'linewidth', 3)
plot(-7.9, -11.85, 'rx', 'markersize', 12, 'linewidth', 3)
axis([-10 -6 -20 20]), hold off
While this method can make the root loci pass through the desired closed-loop poles, the compensator will consequently be improper. Note that the closed-loop poles of the compensated system have six additional poles.
  3 Comments
Sam Chak
Sam Chak on 28 Nov 2023
Yup, @Arpad. The Angle Deficiency Compensation method is commonly taught in the Root Locus to determine the locations of the compensator's zero(s) and pole(s). As shown in yours and @Paul's calculations, introducing the compensator's pole at should make the root locus pass through the desired closed-loop poles (see Figure 1). You will also need to make a gain adjustment to move the closed-loop poles to the desired locations (). However, this method only approximates the desired 2nd-order system (see Figure 2), due to the presence of nearby closed-loop poles. But I think the result is acceptable to your Professor.
Previously, I thought your professor wanted you to drive the closed-loop system to behave "exactly" like the reduced 2nd-order system with two poles at , producing the desired overshoot and settling time. That's why I introduced the Dominant Closed-loop Poles concept in which you can freely assign the 3rd desired closed-loop negative real pole that is far far away from the dominant closed-loop poles ().
figure(1)
figure(2)
Sam Chak
Sam Chak on 28 Nov 2023
By the way, these techniques of placing the closed-loop poles in the desired locations fall under the Pole Placement method. I also forgot to mention that the solution to your problem is not unique because there are infinitely many solutions. You can present alternative solutions to your Professor and ask for opinions (considering the pros and cons of each technique).

Sign in to comment.

Products


Release

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!