deg
|
|
|
|
|
|
|
|
6
|
2
|
0.1
|
24.4
|
3
|
60
|
|
|
2.
Calculations
I. Analyze of uncompensated system
) Calculation of the transfer
function of the whole uncompensated system:) without Matlab
Theoretical information:
A transfer function of the analog system
is the ratio of the Laplace transform of the output signal to the Laplace
transform of the input one under zero initial conditions. A continuous time
SISO (single-input-single-output) transfer function is characterized by its numerator and denominator , both polynomials of the Laplace variable
.
Calculations:
The control system of nuclear reactor rods
position looks like:
. 2
So, the transfer function of open-loop
uncompensated system is:
And the transfer function of closed-loop uncompensated
system is:
b) with Matlab
Program code:
disp ('Task 1. Analyze of uncompensated system')
=tf([5]);=tf([6], [1 6
0]);=tf([0.2]);=tf([10.5]);=tf([1], [0.1 1]);
=series (W1, W2);=series (W12, W3);
('Transfer function of open-loop uncompensated
system is')=series (W123, W4)('Transfer function of closed loop uncompensated
system is')=feedback (Wolun, W5, - 1)% task 1.1
Results of the program:
Task 1. Analyze of uncompensated systemfunction
of open-loop uncompensated system isfunction:
^2 + 6 sfunction of closed loop uncompensated
system isfunction:
.3 s + 63
.1 s^3 + 1.6 s^2 + 6 s + 63
2) define zeroes and poles of
obtained transfer function
Theoretical information:
Poles of the transfer function are the roots of
the system characteristic equation. Characteristic equation is the denominator
of transfer function reduced to zero.of the transfer function are the roots of
equation which is obtained by reducing the numerator of the transfer function
to zero.method is based on using of the operators pole(sys) which
calculates the poles of the transfer function and zero(sys) which
calculates zeros.method is based on using of the operator pzmap(sys), which
plots the pole-zero map of the continuous-time LTI model sys on a
complex plane. For SISO systems, pzmap plots poles and zeros of the
transfer function.
Program code:
Results of the program:
Poles and zeros of closed loop uncompensated
system=
.8199
.5901 + 6.4933i
.5901 - 6.4933i=
.0000
.3
3) plot the time responses for obtained
transfer function; the responses should be produced in two different windows of
one figure graphic object using command subplot;
Theoretical information:
There are two main time responses, namely: a step
response and an impulse one.response of the dynamic system on the step
(Heavyside) function at its input is called a step response. To plot the step
response of an arbitrary LTI (Linear Time Invariant) model sys on the
screen, operator step(sys) is used.
The response of the dynamic system on the -function at its input is called an
impulse response. To plot the impulse response of an arbitrary LTI model sys
on the screen, operator impulse(sys) is used. The duration of simulation
is determined automatically to display the transient behavior of the response.
Program code:
figure(2)(2,1,1), step(Wclun),on, legend ('Step
of uncompensated system')% task 1.3(2,1,2), impulse(Wclun),on, legend ('Impulse
of uncompensated system')
Results of the program:
Fig.4
4) determine settling time and overshoot
of the uncompensated system;
Theoretical information:
Settling time is the period of time from the beginning of transients to the
moment of time, after which the inequality
takes place, where is given small constant value, which is
usually equal to of the steady state value. Or in other
words, it is the time for the system output to settle down to within a
tolerance band of the final value, normally between ±2 or 5%.
Overshoot is maximal deflection of the output value
in the transient process of the transient process from the steady state value.
It normally expressed as percentage determined as:
,
where is maximal deflection of the transient process.the most cases a
system would has sufficient stability margins if overshoot is less or equal to .
Results of the program:
From the Fig.4:of closed loop compensated
system is 81.8%.time of closed loop uncompensated system is 6.76 sec.
5) check whether the system is stable or
not using Hurwitz criterion;
Theoretical information:
Hurwitz stability criterion
The characteristic equation of order system is
.
For this system to be stable it is
necessary and sufficient that the determinant of Hurwitz matrix and the determinants
of all its diagonal minors are positive.matrix is
.
The main diagonal of Hurwitz matrix
contains coefficients of characteristic equation beginning with to . Elements which are located above the main diagonal have
increasing indices; elements which are located under the main diagonal have
decreasing indices. So odd-numbered lines contain coefficients of
characteristic equation with odd indexes, and even-numbered lines contain
coefficients of characteristic equation with even indexes. Places, in which
coefficients are absent, are filled by zeros.
Program code:
[n0, d0]=tfdata (Wclun, 'v')d0 (1)>0&d0
(2)>0&d0 (3)>0&d0 (4)>0&(d0 (2)*d0 (3) - d0 (1)*d0
(4))>0('System is stable by Hyrwitz criterion')('because all coeficients in
Hurwitz matrix are positive')('and particular case for 3rd order s-m is
true')disp ('system is unstable')% task 1.5
Results of the program:
n0 =
0
6.3000 63.0000=
.1000
1.6000 6.0000 63.0000is stable by Hyrwitz criterion because all coeficients in
Hurwitz matrix are positive and particular case for 3rd order s-m is true
6)
check whether the system is stable or not using Nyquist criterion;
Theoretical
information:
There
are two formulations of Nyquist criterion of stability depending on whether the
open-loop system is stable or not stable.the open-loop system is stable, the
closed loop system would be stable if the Nyquist plot of the open-loop system
does not encircle the point with coordinates (-1, 0) at the complex plane.
If the open-loop system is unstable, the
closed loop system would be stable if the Nyquist plot of the open-loop system
encircles the point with coordinates (-1, 0) at the complex plane times, where is the number of roots in RHP, in the
counter-clockwise direction when frequency changes from zero to infinity.
Operator nyquist(sys) plots the Nyquist
response of an arbitrary LTI model sys.
Program code:
figure(3), nyquist(Wolun), grid on% task 1.6
Results of the program:
. 5
System is stable because Nyquist plot for the
open-loop system doesn`t encircle point (-1; 0).
7) if system is stable, determine the
stability margins using Bode diagram;
Theoretical information:
Stability margins
There is a special operator margin for
calculating the stability margins in MatLab. Its syntax is
[Gm, Pm, Wcg, Wcp]=margin(sys),
where is gain margin, is phase margin, is frequency under which , is frequency under which . and are crossover frequencies.invoked without
left-hand arguments, margin(sys) plots the open-loop Bode response with
the gain and phase margins marked by vertical lines.
Program code:
figure(4), margin(Wolun), grid on, legend
('Uncompensated system')% task 1.7
Results of the program:
. 6
So, gain margin is equal to infinity and phase
margin - 41.1 deg.
8) do the conversion from the transfer
functions to state space;
Theoretical information:
State space description
Complete state space description of a continuous
system with constant coefficients is the following:
To create the continuous-time LTI (Linear Time
Invariant) models in MatLab it is necessary to use operator ss.example,
sys=ss (A, B, C, D)
creates the state space of a continuous-time
system.
Relationship between state space and
transfer function description
The diagram explaining the relationship between
the state space description and description with the transfer functions is
represented in fig. 7.
. 7
For the continuous case we have
.
In MatLab operator ss2tf carries out
conversion from state space description to transfer functions. If output is a
vector than numerator Num of the transfer function is a matrix, each row
of which contains numerator polynomial coefficients of the corresponding
transfer function (from input to the first, second and so on outputs). If it is
necessary to define any concrete transfer function then it is necessary to take
the corresponding row from the matrix Num, for example tf2=tf (Num(2,:),
Den). Denominator would be common for all transfer functions.
Program code:
[A0, B0, C0, D0] = tf2ss (n0, d0)
Results of the program:
A0 =
-60 -630
0 0
1 0=
=
63.0000 630.0000=
9) calculate the eigenvalues of the
uncompensated system;
Theoretical information:
Eigenvalues
The eigenvalues of the square matrix are the roots of the equation: , where is a unit matrix. The result of using the operator is a vector containing the eigenvalues of
the square matrix .
Program code:
eigenvalues = eig(Wclun)% task 1.9
Results of the program:
eigenvalues =
.8199
.5901 + 6.4933i
.5901 - 6.4933i
10) calculate -norm of the uncompensated system.
Theoretical information:
Calculation of -norm
The -norm of a stable continuous system with a transfer function is a square root from average value of a
square of impulse response of a system, and at conversion to Laplace
transformation according to Parseval theorem, this norm is determined by the
following way:
,
where is trace of matrix; is matrix conjugated to ; is the integral of the error squared
under input signal which is equal to -function.index can be calculated with the help of the controllability or
observability gramians. The function gram allows to calculate Gram
function for an estimation of system controllability, which is called controllability
gramian. Controllability gramian is applied to research controllability
properties of the system models set in state space, and also for construction
of their minimal realizations. It is more convenient for calculations, than
controllability matrix.gramian of this system is determined with the integral:
.
In MatLab: Gc=gram (A, B)gramian is
positive defined if and only if the pair of matrices is controlled. Controllability gramian is
calculated with the solution of the continuous Lyapunov equation:
.
In MatLab: Gcc=lyap (A, B*B').
Observability gramian of the system is determined
with the integral:
.
In MatLab: Go=gram (A, C).gramian
is positive defined if and only if the pair of matrixes is observed. Observability gramian is
calculated with the solution of the continuous Lyapunov equation:
.
In MatLab: Goo=lyap (A, C’*C).is
that the matrix should be stable, i.e. for continuous
models all eigenvalues should have negative real parts.is illustrated in MATLAB
using the following commands: gram (controllability gramian), lyap
(the solution of Lyapunov equation), trace (a trace of matrix).
Gc=gram (A, B);%or Gc=lyap (A, B*B');
J=sqrt (trace(C*Gc*C'));
%or J=sqrt (trace(B'*Go*B)), where Go is
observability gramian.
-norm can also be determined in MatLab as follow:
h2_n=normh2 (A, B, C, D)
-norm accepts the value which is equal to infinity in the
following cases:
- the model is
unstable or neutral;
- the continuous
model has a nonzero matrix D.
Program code:
H2op_loop = normh2 (A0, B0, C0, D0)% task 1.10
Results of the program:
H2op_loop =
.6137
II. Synthesize the real PD-compensator () which would guarantee desired phase
margin at gain crossover frequency :
Theoretical information:
Synthesis of PD-controllers
Transfer function of ideal PD-controller
is the following:
,
where is the gain of proportional part; is the gain of differential one.gain increased unrestrictedly
results in the infinite gain on high frequencies. So, in order to restrict the
gain on high frequencies an additional pole is appended in the differential
component of PD-controller. In this case the transfer function of real
PD-controller can be written as:
,
where is a very small value (). If is known so it is necessary to determine
two parameters in the synthesis procedure.synthesis of PD-controller, first of
all, it is necessary to calculate the phase shift of controller at frequency as
and to determine
.
.
In MatLab:
theta=-pi+Pm/57.296-phase/57.296=cos(theta)/mag=sin(theta)/(om*mag)
% Creation of the transfer function of
real PD-controller=tf([Kp*tau+Kd Kp], [tau 1])
code:
disp ('Task 2. Synthesis of real PD-compensator')
t1=0.0001;=2.5;=60;
1) calculate the phase shift of compensator at given frequency ;
Program code:
disp ('Phase shift of compensator at given
frequency w ')
[mag, phase] = bode (Wolun, w)% task
2.1=-pi+Pm/57.296-phase/57.296;
Results of the program:
mag =
.8769=
.6199
2) determine the gain of proportional part
of compensator ;
Program code:
disp ('Coefficient of proportional unit
')=cos(theta)/mag% task 2.2
Results of the program:
Coefficient of proportional unit=
.2558
3) determine the gain of differential part
of compensator ;
Program code:
disp ('Coefficient of differentiated unit
')=sin(theta)/(w*mag)% task 2.3
Results of the program:
Coefficient of differentiated unit=
.0133
4) plot Bode diagram of the
PD-compensator.
Program code:
disp ('Tf of compensation unit')=tf([Kp*t1+Kd
Kp], [t1 1])% tf of compensation unit
('Tf of compensator')=Wc% tf of
compensator(5)(Wcomp), grid on,('Compensator')% task 2.4
Results of the program:
Tf of compensation unitfunction:
.01323 s + 0.2558
.0001 s + 1of compensatorfunction:
.01323 s + 0.2558
.0001 s + 1
. 7
III. Analysis of compensated system:
) build in Simulink the compensated
system;
Fig. 8
2) determine settling time and
overshoot of the compensated system using transient process received in
Simulink;
Fig. 9
time - 2.4; overshoot - 23.4%.
3) calculate the transfer function of the
whole compensated system;
Program code:
disp ('Tf of open-loop compensated system')=
series (Wolun, Wcomp)('Tf of closed-loop compensated system')=feedback
(Wolcomp, W5, - 1)% task 3.3
Results of the program:
Tf of open-loop compensated systemfunction:
.8334 s + 16.12
.0001 s^3 + 1.001 s^2 + 6 sof closed-loop
compensated systemfunction:
.08334 s^2 + 0.7781 s + 16.12
e-005 s^4 + 0.1002 s^3 + 1.601 s^2 + 5.167 s +
16.12
4) plot the step response for
uncompensated and compensated system; the responses should be produced in one
window of figure graphic object using command hold on;
Program code:
figure(7)(Wclcomp), grid on,on,
step(Wclun),('Closed-loop compensated', 'Closed-loop uncompensated')% task 3.4
Results of the program:
. 10
5) check whether the system is stable or
not using sufficient condition of stability;
Program code:
disp ('Poles of closed-loop compensated
system')(Wclcomp),(8)(Wclcomp), grid on% task 3.5
Results of the program:
Poles of closed-loop copmensated system=
.0e+003 *
.0000
.0015 + 0.0032i
.0015 - 0.0032i
. 11
So the system is stable as all poles are
negative.
6) if system is stable, determine the
stability margins using Bode diagram;
Program code:
figure(9)(Wolcomp), grid on,('Compensated
system')% task 3.6
Results of the program:
Fig. 12
So gain margin is 17.1 dB and phase margin - 60
deg.
7) do the conversion from the transfer
functions to state space;
Program code:
[n1, d1]=tfdata (Wclcomp, 'v'),
[A1, B1, C1, D1]=tf2ss (n1, d1)% task 3.7
Results of the program:
n1 =
0
-0.0833 0.7781 16.1154=
.0000
0.1002 1.6006 5.1666 16.1154=
.0e+006
*
.0100
-0.1601 -0.5167 -1.6115
.0000
0 0 0
0.0000
0 0
0
0.0000 0=
=
.0e+006
*
-0.0083
0.0778 1.6115=
8)
calculate the eigenvalues of the compensated system;
Program
code:
cl_loop_eigen=eig(A1)%
task 3.8
Results
of the program:
cl_loop_eigen
=
.0e+003
*
.0000
.0130
.0015
+ 0.0032i
.0015
- 0.0032i
9) calculate -norm of the compensated system;
Program code:
H2cl_loop=normh2 (A1, B1, C1, D1)% task 3.9
Results of the program:
H2cl_loop =
.4874
Conclusion
In this term paper we estimated properties of
uncompensated dynamic system, synthesized real PD-compensator and as a result
of including it in our system is changed. Also we estimated properties of
compensated system and compared them with corresponding properties of
uncompensated system. So we can say that overshoot is less and settling time of
compensated system is more than in uncompensated system.
References
uncompensated system phase
1. Dorf R.C., Bishop R.H. Modern Control Systems. - Pearson
Education, Inc., Pearson Prentice Hall, 2008. - 1046 p.
2. Бессекерский
В.А., Попов Е.П. Теория систем автоматического управления. - СПб.: Профессия,
2003. - 752 с.