Matlab and Routh criterion for evaluation of K for stability

50 views (last 30 days)
omar khasawneh
omar khasawneh on 13 Apr 2022
Answered: Paul on 13 Apr 2022
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

Accepted Answer

Paul
Paul on 13 Apr 2022
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 = fliplr(coeffs(e,s));
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 reuslt? 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.

More Answers (0)

Products

Community Treasure Hunt

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

Start Hunting!