Need help finding the stability of a closed loop system and storing the values.

I need to loop through the values of Kd (0<Kd<10) and Kp (0<kp<200), for the transfer function. and then created a new closed loop transfer function defined as GCl. Now when looping through this, and testing for the stability of the system in variable B, the new vector x or y is not storing the stable values of Kd and Kp. Basically, is there a way to loop through the values of Kd and Kp for the transfer function GCl and then output the stable values of Kd and Kp? And lastly how can I plot the stable region?

clc
clear
Gred = tf([7.385,476.7],[1,19.39,491.5,3835,0])
Kp = 0:1:200
Kd = 0:1:10
for j=1:10
    for i=1:200
        C = tf([Kd(j),Kp(i)],[0.005,1]);
        num = Gred*C;
        den=1+Gred*C;
        B = isstable(num,den);
        %Check if closed loop map is stable with given pair
        %If it is stable, then store Kp, Kd in their vectors:
        if B == 0
            x = Kd;
            y = Kp;%Will overwrite the values of Kp and Kd to the ones we need 
        end
        i = i+1;
    end
end
% plot(x,y)

Thank you so much in advance and looking forward to any help!

[Paul]

Before going any further can you clarify if the problem is in continuous time or discrete time? I'm asking because Gred and C are formed for continuous time transfer functions, but that use isstable() looks like it is for the numerator and denominator pair of a discrete time transfer function, but then num and den aren't polynomials as would be required for that use case.

Use the feedback() command to form the closed loop system.

Also, Kp has 201 elements and Kd as 11 elements, so the j and i loops need to be adjusted accordingly.

Answer 1

[Paul]

Assuming continuous time:

Gred = tf([7.385,476.7],[1,19.39,491.5,3835,0]);
Kp = 0:1:200;
Kd = 0:1:10;

Preallocate the matrix for the result. Preallocate to NaN so it can be checked later, if desired, to make sure all elements are filled in.

Stable = nan(numel(Kp),numel(Kd));

Loop over the compensator gains

for j=1:numel(Kd)
    for i=1:numel(Kp)
        C = tf([Kd(j),Kp(i)],[0.005,1]);
        Clsys = feedback(C*Gred,1);
        Stable(i,j) = isstable(Clsys);
    end
end

Edit: There are better ways to visualize this than surf().

surf(Kd,Kp,Stable,'EdgeColor','none'),view(2),colorbar
xlabel('Kd')
ylabel('Kp')
title('Yellow = Stable, Blue = Unstable')

The blue sliver for Kp = 0 and Kd >= 1 comes about because in those cases Clsys has pole and a zero at s = 0. But isstable() only only looks at the poles and considers a pole at the origin to be unstable. if that's not the desired behavior, then change the code to

Clsys = minreal(feedback(C*Gred,1));

[Sam Chak]

+1. Good job, @Paul. It's a very nice and informative plot. By the way, can you insert the texts to display the stable and unstable regions?

Also, is the horizontal bar generated by the surf()?

[Paul]

I couldn't get annotations on the plot to work, so I added a title. I guess the colorbar isn't really necessary with the title added.

Yes, surf() generates the blue sliver at Kp = 0, Kd >= 1, as it should.

[Nimit Chovatia]

This is so helpful, thank you so much @Paul! I apologise for the late gratitude! I did have another question, I was trying to get the settling time out for the stable kp and kd values, any idea how I could do that? I tried this but:

clc 
clear
Gred = tf([7.385,476.7],[1,19.39,491.5,3835,0]);
vecKd = [];
vecKp = [];
for Kp = 0:1:200
    for Kd = 0:0.1:10
    C = tf([Kd,Kp],[0.005,1]);
    Gcl = feedback(Gred*C,1);
    B = isstable(Gcl);
    %Check if closed loop map is stable with given pair
    %If it is stable, then store Kp, Kd in their vectors:
        if B == 1
        vecKd(end+1) = Kd;
        vecKp(end+1) = Kp;
        end
    end
end
plot(vecKd,vecKp,'.')

X=pid(vecKp,0,vecKp);
Z=feedback(X*Gred,1);
t=0:0.1:2;
step(Z,t)

But I don't think it is working as it is taking a very long time to load. I also think it would end up plotting an unstable response as I think it is plotting all the Kp and Kd values instead of just the stable ones. I would appreciate any suggestions! Thank you in advance.

[Paul]

Hi Nimit,

Here is the first part of your code

clc 
clear
Gred = tf([7.385,476.7],[1,19.39,491.5,3835,0]);
vecKd = [];
vecKp = [];
for Kp = 0:1:200

I changed the step size of Kd to meet the execution time constraint of Answers.

    for Kd = 0:10
    C = tf([Kd,Kp],[0.005,1]);
    Gcl = feedback(Gred*C,1);
    B = isstable(Gcl);
    %Check if closed loop map is stable with given pair
    %If it is stable, then store Kp, Kd in their vectors:
        if B == 1
        vecKd(end+1) = Kd;
        vecKp(end+1) = Kp;
        end
    end
end
plot(vecKd,vecKp,'.')

Here is a problem, vecKp is being used in both arguments

% X=pid(vecKp,0,vecKp);
X = pid(vecKp,0,vecKd);

But the compensator used to derive the valid values of Kp and Kd had another term in the denominator.

C = X/(0.005*tf('s') + 1);

Closed loop systems

Z = feedback(C*Gred,1);

It takes too long on Answers to run all of them, so just run the first few.

t=0:0.01:2;
step(Z(1,1,1,1:10),t);

Get the step response information

for ii = 1:10
info(ii) = stepinfo(Z(1,1,1,10));
end
info(5) % for example

ans = struct with fields:
         RiseTime: 25.7201
    TransientTime: 51.5012
     SettlingTime: 51.5012
      SettlingMin: 0.9000
      SettlingMax: 0.9985
        Overshoot: 0
       Undershoot: 0
             Peak: 0.9985
         PeakTime: 93.3688