Control of an electromechanical system based on DPT. Desired LFC method and other Matlab tools

A common task when teaching the theory of automatic control is to calculate the correcting device using the desired LFC method. This problem is given to introduce you to the larger world of frequency domain control.
Why the frequency method at all when there is a modal one?
The fact is that in 1978, John Doyle in an article Guaranteed Stability Margins for LQG Regulators it was shown that for LQG controllers there is no guaranteed stability margin, and therefore, depending on the control object, noise and interference in the control and measurement channels, the LQG controller can be arbitrarily sensitive to uncertainty in the model and time delays, which means it can be however unreliable (robust).
In this article we will show several ways to calculate a compensator using frequency methods, in addition to the desired LFC method, in the Matlab package using the Control System Toolbox.

1 Control object (Plant)

The control object (OU), for which we will develop a compensator, is an electromechanical system that includes a DC electric motor Mpower converter TPamplifier A1filter F and engine speed sensor (permanent magnet tachogenerator) G.

System parameters:

\begin{array}{l} L_P =0\ldotp 01\;\mathrm{Henrys}-\mathrm{inductance}\;\mathrm{converter}\\ R_P =0\ldotp 9\;\mathrm{Om}- \mathrm{resistance}\;\mathrm{converter}\\ L_D =0\ldotp 01\;\mathrm{Henrys}-\mathrm{inductance}\;\mathrm{motor}\\ R_D =4\;\mathrm{ Om}-\mathrm{resistance}\;\mathrm{motor}\\ c_e =1\ldotp 2\;V\cdot s-\mathrm{constructive}\;\mathrm{constant}\;\mathrm{DPT}\ ;\mathrm{by}\;\mathrm{voltage}\;\\ c_M =\;1\ldotp 2\;N\cdot \frac{m}{A}-\mathrm{constructive}\;\mathrm{constant }\;\mathrm{DPT}\;\mathrm{by}\;\mathrm{moment}\\ J_D =0\ldotp 04\;\mathrm{kg}\cdot m^2 -\mathrm{given}\; \mathrm{moment}\;\mathrm{inertia}\;\\ k_{\mathrm{TG}} =0\ldotp 9\;V\cdot s-\mathrm{KU}\;\mathrm{tachogenerator}\\ k_P =10-\mathrm{KU}\;\mathrm{power}\;\mathrm{converter}\\ R_1 =10\;\mathrm{kOm}-\mathrm{resistance}\;\mathrm{OS}\; \mathrm{O-Amp}\\ R_2 =R_3 =20\;\mathrm{kOm}-\mathrm{input}\;\mathrm{resistance}\;\mathrm{O-Amp}\;and\;\mathrm{resistance}\ ;R_3 \\ R_Ф =20\;\mathrm{kOm}-\mathrm{resistance}\;\mathrm{filter}\\ C_Ф =5\ldotp 1\;\mathrm{uF}-\mathrm{capacitance}\; \mathrm{filter} \end{array}

Block diagram of the system in general form:

2 Transfer functions and Laplace transform

Let's write down a list of symbolic variables that we will operate with in further work.

syms K_ampl     % КУ силового преобразователя
syms K_tach     % КУ тахогенератора
syms RF CF TF   % сопротивление, емкость и постоянная времени фильтра
syms R1 R2      % сопротивление ОС и входное, операционного усилителя
syms Rrot Lrot  % сопротивление и индуктивность якоря ДПТ
syms Ce         % конструктивная постоянная ДПТ по напряжению
syms Cm         % конструктивная постоянная ДПТ по моменту
syms b          % коэффициент вязкого трения ДПТ
syms J          % приведенный момент инерции
syms Urot
syms M
syms i
eqn2 = J*diff(omega) + b*omega == Cm*i - M; % механическое уравнение ДПТ
eqns = [eqn1; eqn2] % итоговая система уравнений

\left(\begin{array}{c} \mathrm{Urot}\left(t\right)=\mathrm{Lrot}\,\frac{\partial }{\partial t}\;i\left(t\ right)+\mathrm{Ce}\,\omega \left(t\right)+\mathrm{Rrot}\,i\left(t\right)\\ J\,\frac{\partial }{\partial t }\;\omega \left(t\right)+b\,\omega \left(t\right)=\mathrm{Cm}\,i\left(t\right)-M\left(t\right) \end{array}\right)

2.2 Transfer functions of DBT for disturbance and control

To find the solution to our differential equation, i.e. output signal y\left(t\right)it can be represented as the action of some operator on the input signal u\left(t\right). By Laplace transforming the differential equation under zero initial conditions, we reduce the left and right sides to the form:

W\left(s\right)=\frac{Y\left(s\right)}{U\left(s\right)}

Where Y\left(s\right)=L\left\lbrack y\left(t\right)\right\rbrack, U\left(s\right)=L\left\lbrack u\left(t\right)\right\rbrackA W\left(s\right) – the so-called transfer function.

Laplace transform L time function f\left(t\right)(original) into a function of a complex variable F\left(s\right) (image) is determined using the relation:

L\left\lbrack f\left(t\right)\right\rbrack =\int_0^{\infty } f\left(t\right)\cdot e^{-s\cdot t} \cdot \mathrm{dt }=F\left(s\right)

Where s – variable with dimension c^{-1}.

One useful property of this transformation for us is that the differentiation transformation has the form:

L\left\lbrack \frac{\mathrm{d}}{\mathrm{d}t}f\left(t\right)\right\rbrack =s\cdot F\left(s\right)-f\left (0\right)

And since we agreed that the initial conditions are equal to zero, then any linear differential equation after the Laplace transform will turn into an ordinary polynomial.

Let us transform the first and second differential equation (DE) respectively:

eqnl1 = laplace(eqn1); % преобразование Лапласа уравнения 1
eqnl2 = laplace(eqn2); % преобразование Лапласа уравнения 2
syms s W I u f % переменная s, W = L[omega
eqnl1 = subs(eqnl1,[laplace(Urot
H_M = simplify(subs(lhs(eqnM)/rhs(eqnM)*W/f,u,0))

-\frac{\mathrm{Rrot}+\mathrm{Lrot}\,s}{\mathrm{Rrot}\,b+\mathrm{Ce}\,\mathrm{Cm}+J\,\mathrm{Rrot}\ ,s+\mathrm{Lrot}\,b\,s+J\,\mathrm{Lrot}\,s^2 }

2.3 Transfer functions of the entire system for control and disturbance

Let's write down the transfer functions of each element

Transfer function of the power converter: H_P \left(s\right)=K_{\mathrm{ampl}}.

Tachogenerator transfer function: H_{\mathrm{TG}} \left(s\right)=K_{\mathrm{tach}}.

Filter transfer function: H_Ф \left(s\right)=\frac{1}{\mathrm{CF}\,\mathrm{RF}\,p+1}.

Operational amplifier transfer function: H_U \left(s\right)=\frac{\mathrm{R1}}{\mathrm{R2}}.

H_FA = K_ampl; % передаточная функция силового преобразователя
H_G = K_tach; % передаточная функция тахогенератора
H_F = 1/(RF*CF*s+1); % передаточная функция фильтра
H_OA = R1/R2; % передаточная функция операционного усилителя

Let's compose the transfer function of a closed-loop control system:

W1 = H_OA*H_FA*H_U;
W2 = H_G*H_F;
W_sys_u = simplify(expand(W1/(1+W1*W2))) % передаточная функция по управлению всей системы

\frac{\mathrm{Cm}\,K_{\textrm{ampl}} \,R_1 \,{\left(\mathrm{CF}\,\mathrm{RF}\,s+1\right)}}{ \mathrm{Ce}\,\mathrm{Cm}\,R_2 +R_2 \,\mathrm{Rrot}\,b+J\,\mathrm{Lrot}\,R_2 \,s^2 +\mathrm{Cm} \,K_{\textrm{ampl}} \,K_{\textrm{tach}} \,R_1 +J\,R_2 \,\mathrm{Rrot}\,s+\mathrm{Lrot}\,R_2 \,b\ ,s+\mathrm{CF}\,\mathrm{Ce}\,\mathrm{Cm}\,R_2 \,\mathrm{RF}\,s+\mathrm{CF}\,R_2 \,\mathrm{RF}\ ,\mathrm{Rrot}\,b\,s+\mathrm{CF}\,J\,\mathrm{Lrot}\,R_2 \,\mathrm{RF}\,s^3 +\mathrm{CF}\, J\,R_2 \,\mathrm{RF}\,\mathrm{Rrot}\,s^2 +\mathrm{CF}\,\mathrm{Lrot}\,R_2 \,\mathrm{RF}\,b\ ,s^2 }

Let's compose the transfer function of a closed-loop system based on disturbance:

W3 = W1*-W2;
W4 = 1/(1-W3);
W_sys_M = simplify(expand(H_M*W4)) % передаточная функция по возмущению всей системы

-\frac{R_2 \,{\left(\mathrm{CF}\,\mathrm{RF}\,s+1\right)}\,{\left(\mathrm{Rrot}+\mathrm{Lrot}\ ,s\right)}}{\mathrm{Ce}\,\mathrm{Cm}\,R_2 +R_2 \,\mathrm{Rrot}\,b+J\,\mathrm{Lrot}\,R_2 \,s ^2 +\mathrm{Cm}\,K_{\textrm{ampl}} \,K_{\textrm{tach}} \,R_1 +J\,R_2 \,\mathrm{Rrot}\,s+\mathrm{Lrot }\,R_2 \,b\,s+\mathrm{CF}\,\mathrm{Ce}\,\mathrm{Cm}\,R_2 \,\mathrm{RF}\,s+\mathrm{CF}\,R_2 \,\mathrm{RF}\,\mathrm{Rrot}\,b\,s+\mathrm{CF}\,J\,\mathrm{Lrot}\,R_2 \,\mathrm{RF}\,s^3 +\mathrm{CF}\,J\,R_2 \,\mathrm{RF}\,\mathrm{Rrot}\,s^2 +\mathrm{CF}\,\mathrm{Lrot}\,R_2 \,\ mathrm{RF}\,b\,s^2 }

2.4 Solving the resulting system using the inverse Laplace transform:

Inverse Laplace transform L^{-1} functions F\left(s\right) into a time function f\left(t\right) carried out according to the formula:

f\left(t\right)=L^{-1} \left\lbrack F\left(s\right)\right\rbrack =\frac{1}{2\pi j}\int_{\sigma -j \cdot \infty }^{\sigma +j\cdot \infty } F\left(s\right)\cdot e^{s\cdot t} \cdot \mathrm{ds}

Where j=\sqrt{-1}and constant \sigma is chosen so that the integral converges.

The inverse Laplace transform allows you to conveniently find a solution to the DE under standard input influences for which the originals and images are known.

Let's find a solution to our system when controlling a step signal.

param % загрузка численных параметров системы
syms t % время
g0 = heaviside(t-1); % входной сигнал
G0 = laplace(g0,s); % изображение для входного сигнала
Yp = subs(G0*W_sys_u); % изображение выходного сигнала
ytu = ilaplace(Yp); % оригинал выходного сигнала
fplot(ytu,[0 10])
title('Реакция на ступенчатое управляющее воздействие')
grid on

Let us find the response of the system to a constant disturbing single influence.

g0 = 1; % входной сигнал
G0 = laplace(g0,s); % изображение входного сигнала
Yp = subs(G0*-W_sys_M); % изображение выходного сигнала
ytM = ilaplace(Yp); % оригинал выходного сигнала
fplot(ytM,[0 10])
title('Реакция на единичное возмущающее воздействие')
grid on

Let us find a solution to the system under joint influence of both control and disturbance.

yt = ytu-ytM; % сумма выходных сигналов
fplot(yt,[0 10])
title('Сумма воздействий')
grid on

2.5 Building the system as a Control System Toolbox object

For convenient work with linear systems, the Matlab package provides the Control System Toolbox application.

Let's assemble our system in terms of this application.

Let us find a system of remote control for DBT in the state space.

[V,S] = odeToVectorField(eqns) % уравнения ДПТ в пространстве состояний

\left(\begin{array}{c} -\frac{M\left(t\right)-\mathrm{Cm}\,Y_2 +b\,Y_1 }{J}\\ -\frac{\mathrm{ Ce}\,Y_1 -\mathrm{Urot}\left(t\right)+\mathrm{Rrot}\,Y_2 }{\mathrm{Lrot}} \end{array}\right)

States:

\left(\begin{array}{c} \omega \\ i \end{array}\right)

Vc = feval(symengine,'evalAt',V,'Y = [y1,y2]'); % техническое переименование переменных
syms y1 y2 % символьные переменные
y = [y1;y2];
A = jacobian(Vc,y) % найдем матрицу А модели ДПТ

Matrix A will look like:

\left(\begin{array}{cc} -\frac{b}{J} & \frac{\mathrm{Cm}}{J}\\ -\frac{\mathrm{Ce}}{\mathrm{Lrot }} & -\frac{\mathrm{Rrot}}{\mathrm{Lrot}} \end{array}\right)

A = double(subs(A)); % подставим численные значения в A
B = jacobian(Vc,[Urot;M]) % найдем матрицу B ДПТ

Matrix B:

\left(\begin{array}{cc} 0 & -\frac{1}{J}\\ \frac{1}{\mathrm{Lrot}} & 0 \end{array}\right)

B = double(subs(B)); % подставим численные значения в B
C = [1 0]; % матрица C ДПТ
D = 0; % матрица D ДПТ
DC_Motor = ss(A,B,C,D,'StateName',{'omega','i'},'InputName',{'U','M'},'OutputName',{'omega'}); % модель ДПТ
DC_Motor = zpk(DC_Motor); % передаточная функция ДПТ
DC_Motor.DisplayFormat="time constraint" % формат отображения модели ДПТ

DPT transfer function as a Control System Toolbox object:

DC_Motor =
 
  From input "U" to output "omega":
           0.62176
  --------------------------
  (1+0.09838s) (1+0.004213s)
 
  From input "M" to output "omega":
    -2.5389 (1+0.004082s)
  --------------------------
  (1+0.09838s) (1+0.004213s)
 
Continuous-time zero/pole/gain model.

Let us determine the transfer functions for all elements of the electromechanical system.

Power converter H_P \left(s\right):

H_FA = tf(K_ampl); % передаточная функция силового преобразователя
H_FA.InputName="b"; % вход преобразователя
H_FA.OutputName="U"; % выход преобразователя

Tachogenerator H_{\mathrm{TG}} \left(s\right):

H_G = tf(K_tach); % передаточная функция тахогенератора
H_G.InputName="omega"; % вход тахогенератора
H_G.OutputName="f"; % выход тахогенератора

Filter H_Ф \left(s\right):

H_F = tf(1,[TF 1]); % передаточная функция фильтра
H_F.InputName="f"; % вход фильтра
H_F.OutputName="e"; % выход фильтра

Operational amplifier H_U \left(s\right) :

H_OA = tf(R1/R2); % передаточная функция ОУ
H_OA.InputName="a"; % вход ОУ
H_OA.OutputName="b"; % выход ОУ

Let's assemble a closed system:

sum = sumblk('a = u - e'); % сумматор
Plant = zpk(connect(DC_Motor,H_FA,H_G,H_F,H_OA,sum,{'u','M'},'omega')); % передаточная функция системы
Plant.DisplayFormat="time constraint"; % формат отображения

Let us determine the transfer function of the system based on the control action:

Plant(1)
ans =
 
  From input "u" to output "omega":
                   1.0728 (1+0.102s)
  ---------------------------------------------------
  (1+0.004017s) (1 + 0.1597(0.01906s) + (0.01906s)^2)
 
Continuous-time zero/pole/gain model.

Let us determine the transfer function of the system based on the disturbing influence:

Plant(2)
ans =
 
  From input "M" to output "omega":
           -0.08761 (1+0.004082s) (1+0.102s)
  ---------------------------------------------------
  (1+0.004017s) (1 + 0.1597(0.01906s) + (0.01906s)^2)
 
Continuous-time zero/pole/gain model.

3 LAFCHH and stability margins

Note that when a harmonic signal is applied to the input of a linear system, the output of the system will be the same harmonic signal of the same frequency, but of a different amplitude Aand out of phase \phi. Therefore, to analyze such systems, it is convenient to use the logarithmic amplitude-phase frequency response (LAFC), which demonstrates the system’s response to harmonic signals of different frequencies on a logarithmic scale.

Making a variable change in the transfer function s\to j\cdot w (Where w– frequency), the dependence of the amplitude and phase of the output signal on the frequency of the input signal is determined by the formulas: \left|W\left(j\cdot w\right)\right|=A\left(w\right) And \angle W\left(j\cdot w\right)=\phi \left(w\right).

We can evaluate the quality of a closed-loop system by the frequency characteristics of an open-loop system.

Important assessments of the quality of the system are the stability margins in amplitude and phase. They show how much the system parameters can change in comparison with the calculated ones in order for the system to remain stable.

The amplitude stability margin shows the additional loop gain that is necessary to bring the system to the stability limit.

The phase stability margin shows what additional phase shift is necessary in order to bring the system to the stability limit.

Let us determine the stability margins of our open-loop system

Let's compose the transfer function of an open-loop control system.

Wopen = zpk(connect(H_OA,H_FA,DC_Motor(1),H_G,H_F,'a','e'));
Wopen.DisplayFormat="time constraint"
Wopen =
 
  From input "a" to output "e":
                 27.979
  -------------------------------------
  (1+0.102s) (1+0.09838s) (1+0.004213s)
 
Continuous-time zero/pole/gain model.

Let us estimate the stability margins according to LAFCH.

margin(Wopen)

As can be seen from the graph, the amplitude stability margin Gm = 5.32 dB at a frequency of 69 rad/s and phase stability margin Pm = 9.87 degrees at a frequency of 51 rad/s.

Along with LAFCH, the Nyquist diagram is used to assess stability margins. This is a curve on the complex plane corresponding to the dependence of the real and imaginary parts of the transfer function on frequency. The Nyquist criterion clearly shows the stability of a closed system. If the hodograph of an open-loop system does not cover the point \left(-1;\;0\right)then the system is stable.

Let's construct a Nyquist hodograph of an open-loop system:

nyquist(Wopen)
grid on
xlim([-3 0])
ylim([-1 1])

It is clear from the graph that the hodograph does not cover the point \left(-1;\;0\right)which means the closed system is stable.

The stability margins were found above independently of each other. But in real systems it is often necessary to know what margin there is for simultaneous amplification and delay.
To solve this problem, there is the concept of a marginality disk. This is a circle on the Nyquist diagram, which demonstrates the margin by which the hodograph can still be shifted.
Let's find the permissible simultaneous phase and gain margins:

DM = diskmargin(Wopen); % запасы одновременной устойчивости
figure
h = nyquistplot(Wopen); % диаграмма Найквиста
grid on, hold on
zoomcp(h)
diskmarginplot(DM.GainMargin,'nyquist')
hold off
xlim([-3 0])
ylim([-1 1])

The graph shows that the simultaneous gain and phase margin for the “worst” case is 1.4 dB and 9.5 degrees.

4 Calculation of the corrective link using the desired LFC method

Synthesis of ACS by the desired LFC method is often used in
various student works. This method developed by V.V. Solodovnikov,
is based on the correspondence between logarithmic frequency characteristics
open-loop system and its static and dynamic properties in a closed-loop
condition. The method is used for minimum-phase type systems [Теория автоматического управления под реакцией Воронова А.А. Часть 1, стр. 271, 1986г].

Let us set the following requirements for a closed-loop self-propelled control system:

  • error in steady state \delta =10%percent;

  • transition time t_p =0\ldotp 8\;s;

  • overshoot \sigma =20%percent.

delta = 10; % ошибка
tp = 0.8; % время переходного процесса
sigma = 20; % перерегулирование

Let us construct the asymptotic LFC of the uncorrected open-loop automatic control system.

[L0,~,lgw1,~] = freqasymp(Wopen,wlim,1); % ЛАЧХ ОУ
figure
wlim = [-1 4]; % диапазон интересующих частот
plotL(Wopen,wlim)
hold on

4.1 Desired LFC

Let us construct the desired LFC of an open-loop system that has the desired static and dynamic properties.

The desired LFC consists of three main asymptotes: low-frequency, mid-frequency and high-frequency.

The low-frequency part of the desired LFC determines the static properties. Let's draw this straight line at an angle -20 dB/dec in the low-frequency region of the LFC.

The mid-frequency part determines the stability and margin of stability of the system. Let's carry out the asymptote in the same way with the slope -20 dB/dec crossing the x-axis at the cutoff frequency. The ends of the segment are the place of intersection with the straight lines of the gain margin \Delta L.

Using Solodovnikov’s nomograms we find the cutoff frequency and gain margin.

[omegacut,Lm,~] = nomosol(sigma,tp) % частота среза и запас по усилению соответственно
omegacut = 9.4248
Lm = 28.3741

The high-frequency part is built on the principle of minimal implementation complexity. To do this, we draw the asymptote parallel to the LACCH of the op-amp, i.e. at an angle -60 dB/dec.

As a result, we obtain the following LFC:

4.2 LFC of the correcting device (regulator)

Let's construct the LFC of the sequential correcting device by graphically subtracting the LFC of the op-amp from the desired LFC.

lgw = sort(vertcat(lgwdes,lgw1)); % частоты желаемой ЛАЧХ и ЛАЧХ ОУ
Ldes = interp1(lgwdes,Ldes,lgw); % амплитуды желаемой ЛАЧХ
L0 = interp1(lgw1,L0,lgw); % амплитуды ЛАЧХ ОУ
Lreg = Ldes-L0; % вычитаем из желаемой ЛАЧХ ЛАЧХ ОУ
plot(lgw,Lreg,'LineWidth',2)
hold off
legend('ЛАЧХ ОУ','','','','желаемая ЛАЧХ','+Lm','-Lm','','ЛАЧХ регулятора')

4.3 Determination of the transfer function of the correcting link

Let's calculate the desired gain of the open-loop system.

K0 = 1/(delta/100); % желаемый КУ

Let us determine the transfer function of the correcting link by the type of LFC. A change of -20 dB/dec corresponds to the multiplication of the transfer function by the aperiodic link, and a change of +20 dB/dec corresponds to the multiplication by the inverse aperiodic link. We find the gain depending on the required accuracy.

w1 = 10; % частота первого динамического звена
w2 = 10^2.4; % частота второго динамического звена
T1 = 1/w1; % постоянная времени 1 звена
T2 = 1/w2; % постоянная времени 2 звена
s = tf('s');
KP = K_ampl*K_tach*R1/R2*0.62176; % КУ разомкнутой системы
Kdes = K0/KP; % КУ компенсатора
Compensator = Kdes/s*(T1*s+1)^2/(T2*s+1); % передаточная функция компенсатора
Compensator.InputName="cI";
Compensator.OutputName="cO";
Compensator = zpk(Compensator);
Compensator.DisplayFormat="time constraint"
Compensator =
 
  From input "cI" to output "cO":
  0.35741 (1+0.1s)^2
  ------------------
   s (1+0.003981s)
 
Continuous-time zero/pole/gain model.

4.4 Check

Let's construct the transition function of a closed-loop automatic control system with a compensator.

H_OA2 = H_OA;
H_OA2.InputName="cO";
sum2 = sumblk('cI = u - e');
Sys1 = tf(connect(DC_Motor,Compensator,H_FA,H_G,H_F,H_OA2,sum2,{'u','M'},'omega')); % замкнутая САУ
step(Sys1(1))

Based on the figure, we will determine the main indicators of the quality of regulation. Error in steady state \delta =11% percent, which is slightly more than the required value, but it can be reduced by adjusting the additional amplifier at the input. Transition time t_p =0\ldotp 174\;s. Overshoot \sigma =8% percent.

We will find the gain and phase margins below:

5 Automatic calculation of compensator in Control System Designer

Control System Toolbox in the Matlab package provides us with a convenient application with an intuitive interface and many useful tools for developing controllers – Control System Designer.

To start working in it, you need to prepare transfer functions for an open-loop system without OS and a transfer function of OS.

G = connect(H_OA,H_FA,DC_Motor(1),'a','omega'); % передаточная функция петли
H = connect(H_G,H_F,'omega','e'); % передаточная функция обратной связи

Launch Control System Designer from the APPS panel and select the prepared transfer functions in the Edit Architecture window.

Next, we select a method for calculating the compensator based on the given desired shape of the open-loop system. This method calculates the controller so that the final LFC is as close as possible to the specified one. You can set the LFC using the frequency bandwidth \omega_bthen the LFC will tend to the LFC of the integrator \frac{\omega_b }{s}or LFC of a user-specified transfer function.

Let's compare the quality of regulation of the automatically calculated compensator with the first compensator, which we built using the desired LFC method.

load('Comensators.mat') % загрузим LTI модель компенсатора, подготовленного в Control System Designer
CompensatorCD.InputName="cI";
CompensatorCD.OutputName="cO";
CompensatorCD.DisplayFormat="time constraint"
CompensatorCD =
 
  From input "cI" to output "cO":
  0.32107 (1 + 1.934(0.1071s) + (0.1071s)^2)
  ------------------------------------------
                s (1+0.0175s)
 
Name: C
Continuous-time zero/pole/gain model.

Based on the figure, we will determine the main indicators of the quality of regulation. Error in steady state as well \delta =11% percent, transition time t_p =0\ldotp 153\;sovershoot \sigma =10% percent.

Let's find the gain and phase margins:

6 Automatic calculation of the PID controller in Matlab

Our electromechanical system can also be controlled using a PID controller.

Control System Toolbox has a wide range of
tools for calculating PID controllers. One of the ways is team calculation
pidtune.

pidComp0 = pid(1,1,1); % создадим модель ПИД контроллера
opts = pidtuneOptions('PhaseMargin',74,'DesignFocus','reference-tracking'); % зададим параметры для расчета
[PIDComp,~] = pidtune(Wopen,pidComp0,6,opts) % расчитаем ПИД регулятор для частоты среза 6 рад/с
PIDComp =
 
             1          
  Kp + Ki * --- + Kd * s
             s          

  with Kp = 0.0359, Ki = 0.245, Kd = 0.00131
 
Continuous-time PID controller in parallel form.

Let's compare the transition functions:

Error \delta =11% percent, transition time t_p =0\ldotp 58\;sovershoot \sigma =4% percent.

Gain and phase margins:

7 Automatic calculation of the compensator in the Control System Tuner

Another convenient way to find a compensator that satisfies your parameters is to use the systune command. To do this, you need to prepare a genss model in which the desired compensator will be located.

CompensatorCT = tunableTF('C1',2,2); % искомый компенсатор с 2 нулями и 2 полюсами
CompensatorCT.InputName="cI";
CompensatorCT.OutputName="cO";
Sys4 = connect(DC_Motor,CompensatorCT,H_FA,H_G,H_F,H_OA2,sum2,{'u','M'},'omega'); % genss модель

Next, we will set the target transition function for control and the system’s response to disturbance.

Strack = TuningGoal.StepTracking('u','omega',0.1,10); % зададим требуемую переходную функцию
Srej = TuningGoal.StepRejection('M','omega',1,1); % зададим ограничения по возмущению

Let's do the calculation.

options = systuneOptions('Display','off','MinDecay',1e-9); % выбор опций
[SysTune,~,~] = systune(Sys4,Strack,Srej,options); % поиск удовлетворяющего требованиям компенсатора

Calculation result.

From the figure we determine that the error is in steady state \delta =11\;%percent, transition time t_p =0\ldotp 5\;sovershoot \sigma =0% percent.

Found compensator.

CompensatorCT =
 
  1374.1 (1 + 1.816(0.09973s) + (0.09973s)^2)
  -------------------------------------------
           (1+1.882e-06s) (1+8436s)
 
Name: C1
Continuous-time zero/pole/gain model.

Gain and phase margins:

8 Comparison of expansion joints

Let's compare the operation of an electromechanical system with all 4 compensators with each other with a stepped input signal and a single disturbance.

tend = 10; % время симуляции
st = 0.001; % шаг интегрирования
t = (0:st:tend)';
u = @step_u; % входной сигнал
u_sim = zeros(size
for i = 1:length
    u_sim(i) = u(t(i));
end
Msim = ones(size
figure
lsim(Plant,Sys1,Sys2,Sys3,SysTune,[u_sim Msim],t) % симуляция моделей
legend
xlim([0 3])
ylim([-2 2])

Note that the last 4 compensator turned out to be the “smoothest”, and the first one was the “fastest”. You can see how much better the work with the compensator has become compared to the original system.

The article did not demonstrate all the capabilities of Matlab for calculating regulators. But these examples can significantly make life easier for both students in specialties related to TAU and ordinary engineers.

Link to sources

Link to lectures by Steve Brunton, on which the article was partly written

Similar Posts

Leave a Reply

Your email address will not be published. Required fields are marked *