-
Notifications
You must be signed in to change notification settings - Fork 0
/
rhStabilityCriterion_jpn.m
125 lines (103 loc) · 3.66 KB
/
rhStabilityCriterion_jpn.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
%% Routh-Hurwitz stability criterion
%
% The Routh-Hurwitz stability criterion is a necessary (and frequently
% sufficient) method to establish the stability of a single-input,
% single-output(SISO), linear time invariant (LTI) control system.
% More generally, given a polynomial, some calculations using only the
% coefficients of that polynomial can lead us to the conclusion that it
% is not stable.
% Instructions
% ------------
%
% in this program you must give your system coefficients and the
% Routh-Hurwitz table would be shown
%
% Farzad Sagharchi ,Iran
% 2007/11/12
%% Initialization
% clear ; close all; clc % by 川田
% Taking coefficients vector and organizing the first two rows
% coeffVector = input('input vector of your system coefficients: \n i.e. [an an-1 an-2 ... a0] = '); % by 川田
fprintf('\n')
fprintf('---------------------\n')
fprintf('特性多項式の係数:\n')
coeffVector = denP % by 川田
ceoffLength = length(coeffVector);
rhTableColumn = round(ceoffLength/2);
% Initialize Routh-Hurwitz table with empty zero array
rhTable = zeros(ceoffLength,rhTableColumn);
% Compute first row of the table
rhTable(1,:) = coeffVector(1,1:2:ceoffLength);
% Check if length of coefficients vector is even or odd
if (rem(ceoffLength,2) ~= 0)
% if odd, second row of table will be
rhTable(2,1:rhTableColumn - 1) = coeffVector(1,2:2:ceoffLength);
else
% if even, second row of table will be
rhTable(2,:) = coeffVector(1,2:2:ceoffLength);
end
%% Calculate Routh-Hurwitz table's rows
% Set epss as a small value
epss = 0.01;
% Calculate other elements of the table
for i = 3:ceoffLength
% special case: row of all zeros
if rhTable(i-1,:) == 0
order = (ceoffLength - i);
cnt1 = 0;
cnt2 = 1;
for j = 1:rhTableColumn - 1
rhTable(i-1,j) = (order - cnt1) * rhTable(i-2,cnt2);
cnt2 = cnt2 + 1;
cnt1 = cnt1 + 2;
end
end
for j = 1:rhTableColumn - 1
% first element of upper row
firstElemUpperRow = rhTable(i-1,1);
% compute each element of the table
rhTable(i,j) = ((rhTable(i-1,1) * rhTable(i-2,j+1)) - ....
(rhTable(i-2,1) * rhTable(i-1,j+1))) / firstElemUpperRow;
end
% special case: zero in the first column
if rhTable(i,1) == 0
rhTable(i,1) = epss;
end
end
%% Compute number of right hand side poles(unstable poles)
% Initialize unstable poles with zero
unstablePoles = 0;
% Check change in signs
for i = 1:ceoffLength - 1
if sign(rhTable(i,1)) * sign(rhTable(i+1,1)) == -1
unstablePoles = unstablePoles + 1;
end
end
% Print calculated data on screen % by 川田
% fprintf('\n Routh-Hurwitz Table:\n')
fprintf('\n')
fprintf('---------------------\n')
fprintf('ラウス表:\n')
rhTable
% Print the stability result on screen
fprintf('\n')
fprintf('---------------------\n')
fprintf('安定性の判別:\n')
if unstablePoles == 0
% fprintf('~~~~~> it is a stable system! <~~~~~\n') % by 川田
fprintf(' 安定である\n');
else
% fprintf('~~~~~> it is an unstable system! <~~~~~\n') % by 川田
fprintf(' 安定ではない\n');
end
% fprintf('\n Number of right hand side poles =%2.0f\n',unstablePoles)
fprintf('\n')
fprintf('---------------------\n')
fprintf('不安定極の個数:\n')
fprintf('%2.0f\n',unstablePoles) % by 川田
% reply = input('Do you want roots of system be shown? Y/N ', 's'); % by 川田
% if reply == 'y' || reply == 'Y'
% sysRoots = roots(coeffVector);
% fprintf('\n Given polynomial coefficients roots :\n')
% sysRoots
% end