How to solve Routh Hurwitz with constant K

I am attempting to solve for system stab1ility for the equation: s^4+19*s^3+111*s^2+189*s+K*s+5*K=0

The constant K is throwing me off. I have an idea how to solve this with one variable "s" but need help on how to insert the K as a constant in matlab. thank you in advanced

[Geoff Hayes]

Justin - is K a constant that you are trying to solve for? How would you solve this problem if you didn't have the K?

[justin nabbs]

I normally use the characteristic equation 1+G(s)*H(s)=0 to solve for the equation s^4+19*s^3+111*s^2+189*s+K*s+5*K=0

I then put the constants in the Routh Hurwitz formula to solve. Normally the constants i'm working with are just numbers. I don't know how to define K in matlab so that I can put K in the Routh Hurwitz formula

[justin nabbs]

Solving for stability using Routh Hurwitz gives you the b1,b2 etc. But how do i enter the constant K when i'm entering the coefficients of a characteristic equation through Routh Hurwitz because it gives me a K is undefined error. My equation is s^4+19*s^3+111*s^2+189*s+ K *s+5*K=0 and I used the following syntax:

%%routh hurwitz criteria
clear
clc
%%firstly it is required to get first two row of routh matrix
e=input('enter the coefficients of characteristic equation: ');
disp('-------------------------------------------------------------------------')
l=length(e);
m=mod(l,2);
if m==0
    a=rand(1,(l/2));
    b=rand(1,(l/2));
    for i=1:(l/2)
        a(i)=e((2*i)-1);
        b(i)=e(2*i);
    end
else
    e1=[e 0];
    a=rand(1,((l+1)/2));
    b=[rand(1,((l-1)/2)),0];
    for i=1:((l+1)/2)
        a(i)=e1((2*i)-1);
        b(i)=e1(2*i);
    end
end
%%now we genrate the remaining rows of routh matrix
l1=length(a);
c=zeros(l,l1);
c(1,:)=a;
c(2,:)=b;
for m=3:l
    for n=1:l1-1
        c(m,n)=-(1/c(m-1,1))*det([c((m-2),1) c((m-2),(n+1));c((m-1),1) c((m-1),(n+1))]);
    end
end
disp('the routh matrix:')
disp(c)
%%now we check the stablity of system
if c(:,1)>0
    disp('System is Stable')
else
    disp('System is Unstable');
end

Answer 1

[Julius]

Hi, maybe a bit late, but anyway here is my solution using Matlab and Routh criterion for evaluation of K for stability (root locus does it perfectly in a graphical way by showing critical value of K if locus crosses jw axis or whatever)

syms Kp s
G = (5*s + 2)/(s*(s^2 + 3*s + 2)) % plant TF
Gc = (Kp*(s + 3))/(s + 5) % controller TF
chareq = 1+G*Gc==0
cheq = expand(simplify(chareq))
% haven't figure out how to extract char equation from symbolic, but you can simply copy coefs
% or adapt to you existing char.eq.

%% Update in 2020 - have figured out how to extract coefficients out of char eq
[n, d] = numden(cheq)
cheq = n == 0
cheq = collect(n,s) == 0

R = coeffs(n,s)
%% Now coefficients can be accessed from the vector R and put into Rouht Table
% RT = [R(1,5) R(1,3) R(1,1);
%       R(1,4) R(1,2) 0]

% Routh table first two rows from coefs of char.eq. (from cheq)
RT = [1 17+5*Kp 6*Kp;
      8 10+17*Kp 0];
% the rest of the table
b1 = (RT(2,1)*RT(1,2)-RT(1,1)*RT(2,2))/RT(2,1);
b2 = (RT(2,1)*RT(1,3)-RT(1,1)*RT(2,3))/RT(2,1);
b3 = 0;
c1 = (b1*RT(2,2)-RT(2,1)*b2)/b1;
c2 = (b1*RT(2,3)-RT(2,1)*b3)/b1;
c3 = 0;
d1 = (c1*b2-b1*c2)/c1;
d2 = (c1*b3-b1*c3)/c1;
d3 = 0;
% full Routh table
RT = [1 17+5*Kp 6*Kp;
      8 10+17*Kp 0;
      simplify(b1) b2 b3;
      simplify(c1) c2 c3;
      simplify(d1) d2 d3]
% coeficient Kp values for stability to satisfy condition when b1=0, c1=0 and d1=0
K1 = vpasolve(b1, Kp)
K2 = vpasolve(c1, Kp)
K3 = vpasolve(d1, Kp)

Haven't checked for limitations and what if there is row of zeros or zero in 1st column.

[Dheerendra Kumar]

it's not working

[Rik]

@Dheerendra Kumar how do you expect anyone to solve your problem without information about what is going wrong? If you provide details about what is going wrong (what specific error you get, or how the result is different from what you expect), the chance of getting an answer will greatly increase.

[Muhammad]

This program is not runing, giving an error.

[Julius]

Well.. it would be more informative if you show what kind of error.

[Ahmad]

cheq =

6*Kp + 10*s + 17*Kp*s + 5*Kp*s^2 + 17*s^2 + 8*s^3 + s^4 == 0 & s ~= 0 & s ~= -1 & s ~= -2 & s ~= -5

Error using symengine

Arithmetical expression expected.

Error in sym/mupadmexnout (line 1092)

out = mupadmex(fcn,args{:});

Error in sym/numden (line 15)

[n,d] = mupadmexnout('symobj::numden',s);

Error in My_routh (line 10)

[n, d] = numden(cheq)

[Julius]

remove & s ~= 0 & s ~= -1 & s ~= -2 & s ~= -5 and it should work.

[Ertugrul Bayir]

remove what?

[Julius]

Can't you read above what was the question?

In Ahmad question error "arithmetical expression expected" means that he left strange row of symbols "& s ~= 0 & s ~= -1 & s ~= -2 & s ~= -5" after "==0" and those symboles weren't present in my code.

C'mon.

[fatima mouffok]

chareq = 1+G*Gc==0 remove ==0

chreq = 1+G*Gc like that it's gonna work!!

[ami raa]

thank you so much

Answer 2

[J. Carlos Aguado]

I assume that the origin of this K is a proportional controller. Therefore, the simplest way (according to MatLab implementations, not to mathematical simplicity) to study its effect on stability is to use the root locus of the plant. MatLab will draw it for you and give you the gain value that marks the frontier with instability

[Ahmad Ghareeb]

How to solve it using Routh?