Matlab and Routh criterion for evaluation of K for stability

535 views (last 30 days)
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
Edited: Paul on 4 Jan 2024
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.
  2 Comments
Scott
Scott on 30 Jun 2023
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
Paul on 1 Jul 2023
Hi Scott,
Is there anything materially different between returning parameters and conditions in separate variables as opposed to fields in a struct?

Sign in to comment.

More Answers (0)

Products

Community Treasure Hunt

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

Start Hunting!