Matlab and Routh criterion for evaluation of K for stability

I am traying to evaluation & solve this system stability for the equation s^3 + 3ks^2 + ( k+2 )s + 4 = 0

I strugling with the k constant .

but need help on how to insert the K as a constant for the ROUTH HURWITZ in matlab .

%%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

[Carlos M. Velez S.]

Good job. See function in this link.

Answer 1

[Paul]

Hi @omar khasawneh,

Without a value for K I think you'll need to use the Symbolic Math Toolbox. Your code adapts easily. I did not verify that the code correclty implements the Routh test.

%firstly it is required to get first two row of routh matrix
%e=input('enter the coefficients of characteristic equation: ');
%disp('-------------------------------------------------------------------------')
e = str2sym('s^3 + 3*k*s^2 + (k+2)*s +4')

e = 

syms s
e = coeffs(e,s,'all')

e = 

l=length(e);
m=mod(l,2);
if m==0
    a=sym(rand(1,(l/2)));
    b=sym(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=sym(rand(1,((l+1)/2)));
    b=sym([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=sym(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:')

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
syms k real
ksol = solve(c(:,1)>0,k,'ReturnConditions',true)

ksol = struct with fields:
             k: x
    parameters: x
    conditions: 21^(1/2)/3 - 1 < x

We see that the Routh criteria came up with the solution K > sqrt(21)/3 - 1 = 0.5275

Can we veriy that result? e(s) can be rewritten as:

e(s) = s^3 + 2*s + 4 + k*(3*s^2 + s)

Then we use the allmargin command, which assumes unity feedback, i.e., k = 1

allmargin(tf([3 1 0],[1 0 2 4]))

ans = struct with fields:
     GainMargin: 0.5276
    GMFrequency: 1.5897
    PhaseMargin: [-27.4947 78.6316]
    PMFrequency: [1.1423 3.5586]
    DelayMargin: [5.0802 0.3856]
    DMFrequency: [1.1423 3.5586]
         Stable: 1

The Stable result tells us that with unity feedback (k = 1) the closed loop system is stable, i.e., the roots of e(s) are in the left half plane. The GainMargin result tells us the gain k can be reduce to 0.5276 until instability, which agrees nicely with the symbolic result.

[Scott]

Hi , you were big help for me .

for workiing this code correctly we should use this:

[ksol , parameters , conditions] = solve(c(:,1)>0,k,'ReturnConditions',true) instead of

ksol = solve(c(:,1)>0,k,'ReturnConditions',true)

thank you .

[Paul]

Hi Scott,

Is there anything materially different between returning parameters and conditions in separate variables as opposed to fields in a struct?