Математическая модель объекта управления
Для оценки динамических характеристик сформируем математическую модель
объекта управления.
,
,
Z0=0, Z3=Lz.
Граничные
условия для поверхностей S2, S5 записываются в виде:
Условия
на границах раздела сред, отражающие равенство температур и тепловых потоков,
записываются соотношениями
Управляющее
воздействие в виде теплового потока распределено по границе S1:
Поверхность S6
теплоизолирована
где
- температурное поле i-ой среды; U(x,y,τ)- управляющее воздействие; x,y,z-
пространственные координаты, τ -
время.
Теплофизические параметры заданы следующими значениями (параметры заданы
в системе СИ):
.3 Анализ объекта управления
Рассмотренная выше математическая модель объекта допускает аналитическое
решение (аналитически могут быть определены динамические характеристики объекта
управления по выбранным пространственным модам входного воздействия).
Рассматриваемый объект принадлежит к классу пространственно-инвариантных. В
качестве собственных вектор-функций (пространственных мод) выберем функции вида
Вид собственных вектор функций оператора объекта обусловлен граничными
условиями (2), (3).
Определим реакцию объекта на выбранные моды входного воздействия. Реакция
объекта на выбранную пространственную моду входного воздействия может быть
представлена в виде:
Преобразуя по Лапласу, при нулевых начальных условиях функцию выхода и
входное воздействие и взяв их отношение, получим передаточную функцию
рассматриваемого объекта по выбранной пространственной моде. В рассматриваемом
случае эта передаточная функция может быть записана в виде
.
Записывая передаточную функцию рассматриваемого объекта с использованием
обобщенной координаты, получим
,
В
рассматриваемом случае поставленная задача решалась численно. Для этого,
используя математическую модель объекта, была составлена численная модель и
определена реакция объекта на выбранные пространственные моды входного
воздействия (определена функция для
выбранных значений η
и γ). Схема дискретизации объекта управления приведена на рисунке 2.
Рисунок
2 - Схема дискретизации
При
моделировании объекта управления были выбраны следующие значения переменных:
С(G)
=1000; Nx=5; Ny=5; Nz =15;
Δx=Lx/(Nx-1);
Δy=Ly/(Ny-1);
Δz=Lz/(Nz-1).
Как
известно, в методике синтеза распределенных регуляторов используют динамические
характеристики двух пространственных мод.
В
результате моделирования определена функция и . Их графики приведены на рисунках 3 и 4
соответственно.
Аппроксимируем
передаточную функцию по выбранным пространственным модам передаточной функцией
вида
,
Рисунок
3 - График функции
Рисунок
4 - График функции
В
результате численного моделирования получены следующие значения параметров
передаточной функции:
η=1,γ=1, G1=123.370055, k(G1)=0.076, T(G1)=4076, τ1(G1)=700;
η=3,γ=3 , G3=3084.251375, k(G3)=-0.021 , T(G3)=2753, τз(G3)=403
2. Постановка задачи управления
Постановка задачи: для системы управления объектом, передаточные функции
которого по выбранным пространственным модам заданы в виде (4), синтезировать
распределенный высокоточный регулятор.
При
этом на запасы устойчивости разомкнутой системы и на параметры наложены ограничения: запас устойчивости по фазе ; запас устойчивости по модулю ΔL ≥ 10 дб., значение
параметра =0 .
Передаточная
функция распределенного высокоточного регулятора записывается в виде:
3. Синтез распределенного высокоточного регулятора (РВР)
.1 Методика синтеза РВР
Методика синтеза РВР распадается на следующие этапы:
. Для двух выбранных пространственных мод (G1 и G3) определим желаемые точки среза модуля разомкнутой системы.
При этом положим, что фазовый сдвиг, вносимый в систему регулятором равен нулю.
где
W(G,jω) - комплексный передаточный коэффициент объекта управления.
Используя
уравнение (5), для выбранных пространственных мод (G1 и G3),
определим значения ω1 ω3.
2.
Определение параметров пространственно-усилительного звена
Подставляя
, в
соотношение
определим
значения модуля объекта управления для выбранных пространственных мод. Так как являются частотами среза модуля разомкнутой системы,
то коэффициенты усиления регулятора в этих точках равны:
,
Определение
параметров и будем
осуществлять, исходя из условия
,
,
.
Поделив
(8) на (7), придем к следующему результату:
.
Подставляя
вычисленное значение в (1.23) и преобразуя, получим
.
Определение параметров пространственно - интегрирующего и пространственно -
дифференцирующего звеньев.
Определение
параметров регулятора будем осуществлять, исходя из условия, что значение
частот принадлежит линии перегиба. Для частот, принадлежащих
линии перегиба, фазовый сдвиг, вносимый в разомкнутую систему регулятором,
равен нулю. Подставляя , в
уравнение линии перегиба, получим следующую систему уравнений:
,
.
Вычитая
из (12) (11), придем к следующему результату:
.
Используя
(13), определим значения и .
.1.
Если , то положим .
Тогда
определяется из соотношения
при
этом на изменение значения наложено
ограничение .
Определим
взаимосвязь параметров рассматриваемых звеньев с параметром Δ.
, , ,
Преобразуя,
придем к следующему результату
Подставляя
(16) в (11) и преобразуя, придем к следующему результату
.2.
Если , то положим в (11)-(13) .
Тогда
определяется из соотношения
при
этом на изменение значения наложено
ограничение .
Преобразуя
(16), получим
Подставляя
(19) в (11) и преобразуя, получим соотношение, для вычисления параметра Е2
3.2 Синтез распределенного регулятора для системы управления
температурным полем многослойной пластинки
В соответствии с методикой синтеза РВР, рассмотренной выше:
. Для двух выбранных пространственных мод =123.370055
G3=3084.251375
определим желаемые точки среза модуля разомкнутой системы.
Для определения частот среза модуля разомкнутой системы, в соответствии с
(5), получим следующее уравнение:
, ()
Подставляя
значение , и в (21), определим значение частот среза модуля.
.
Определение параметров пространственно-усилительного звена
Подставляя
, в
соотношение (4), определим значения модуля объекта управления для выбранных
пространственных мод. Так как являются
частотами среза модуля разомкнутой системы, то вычисленные коэффициенты
усиления регулятора в этих точках равны:
,.
Подставляя
вычисленные значения в (9) и (10) получим:
E1= , n1=
.
Определение параметров пространственно - интегрирующего и пространственно -
дифференцирующего звеньев.
Подставляя
вычисленные значения в , получим Δω2>1 .
Подставляя
вычисленные значения Δω2,
G(i) в (14) , (16), (17) получим:
, n4= , E4=, E2= .
Передаточная
функция синтезированного регулятора записывается в виде:
4. Анализ работы замкнутой системы управления
Подавая
на вход РВР входное воздействие , на
входе будем иметь:
;.
Дискретный
аналог алгоритма управления (23) имеет вид:
где
,,
1<ν<Nx , 1
<ξ<Ny .
- точки
дискретизации по оси x и у соответственно; - шаги
дискретизации по осям х и у. Используя передаточную функцию синтезированного
регулятора, запишем алгоритм управления во временной области
1<ν<Nx , 1
<ξ<Ny
,
заданное
значение температурного поля,
текущее
значение температурного поля.
Используя
численные модели объекта и регулятора, осуществлено моделирование работы
замкнутой системы управления. По результатам моделирования построены графики
управляющего воздействия и функции рассогласования в точках ν=3, ξ=3 (см. рисунок 5). При этом
Рисунок
5 - График переходного процесса
Заключение
В настоящей работе в качестве объекта управления выступает температурное
поле многослойной пластинки. В результате выполнения данной курсовой работы мы
изучили динамику его поведения. Нами была построена математическая модель
объекта управления, а также была сформулирована задача управления.
Для построенной модели объекта нами был синтезирован распределенный
высокоточный регулятор, и были построены переходные процессы модели с
включенным в неё регулятором. Синтезированный нами регулятор обеспечивает
плавный нагрев многослойной пластинки до заданной температуры.
Таким образом, поставленная задача выполнена - поведение объекта управления
исследовано, устройство его регулирования рассчитано.
1. А.В. Малков, И.М. Першин - «Системы с распределенными
параметрами. Анализ и синтез».
. Учебное пособие к выполнению курсовой работы по дисциплине
«Моделирование распределенных систем управления».
Приложение
Листинг программы
all= 10;
. Зададим исходные данные:
Геометрические параметры пластины:
Lx = 0.1*NumVar;= 0.2*NumVar;= 0.04*NumVar;
% Слои пластины:= 0.6*Lz;= 0.7*Lz;_zv = 0.8*Lz;
Коэффициент температуропроводности:
a1 = NumVar*1*10^-5;= NumVar*2*10^-4;= 5*a1;
% Дискретизация:= 5;= 5;= 15;=
Lx/(Nx-1);= Ly/(Ny-1);= Lz/(Nz-1);= (6/15)*(Lz/5);= (3/15)*(Lz/2);=
(6/15)*(Lz/5); = 0.01;
Технические параметры:
lambtp1
= NumVar*0.02;
lambtp2 = NumVar*5;=0;=3;= zeros(Nx,Ny);_g = zeros(Nx,Ny);=
zeros(Nx,Ny);= zeros(Nx,Ny,Nz);_etta_gamma_ksi = zeros(Nx,Ny,Nz);
U = zeros(Nx,Ny);
. Расчет данных:
Расчет управления:= 1;= 1509.977;= 1.39;
E2 = 934.1678;= 0.00099;= 80.765158;= 57.231289;=
34.8664826;= 9.278370;= 0.05944;= zeros(5,5);_prev = zeros(5,5);lmb = 1:1*10^5=
t + delt;etta = 2:Nx-1gamma = 2:Ny-1_g(etta,gamma) = 10; etta = 1:Nxgamma =
1:Ny (etta,gamma) = T_g(etta,gamma) - T(etta,gamma,4)etta = 2:Nx-1gamma =
2:Ny-1(etta,gamma) = (F(etta-1,gamma)-2*F(etta,gamma)+F(etta+1,gamma))/delX^2+...
(F(etta,gamma-1)-2*F(etta,gamma)+F(etta,gamma+1))/delY^2;_prev(etta,gamma)
= E4*(((N4-1)/N4)*F(etta,gamma)-(1/N4)*lapl(etta,gamma))*delt;(etta,gamma) =
SUMM(etta,gamma)+SUMM_prev(etta,gamma);
% U(etta,gamma) = 1000*cos(pi*(i-0.5)*etta/(Nx-1))*sin(pi*gamma*(i-0.5)/(Ny-1));(etta,gamma)
=E1*(((N1-1)/N1)*F(etta,gamma)-(1/N1)*lapl(etta,gamma)) + SUMM(etta,gamma) +
E2*((F(etta,gamma)-Fold(etta,gamma))/delt);
end
Тепловые процессы, протекающие в печи
for etta = 2:Nx-1gamma = 2:Ny-1ksi = 2:6_etta_gamma_ksi(etta,gamma,ksi)
= a1*delt*((T(etta-1,gamma,ksi)-2*T(etta,gamma,ksi)+T(etta+1,gamma,ksi))/delX^2
+ ...
(T(etta,gamma-1,ksi)-2*T(etta,gamma,ksi)+T(etta,gamma+1,ksi))/delY^2
+ .
(T(etta,gamma,ksi-1)-2*T(etta,gamma,ksi)+T(etta,gamma,ksi+1))/delZ1^2);etta
= 2:Nx-1gamma = 2:Ny-1ksi = 8_etta_gamma_ksi(etta,gamma,ksi) =
a2*delt*((T(etta-1,gamma,ksi)-2*T(etta,gamma,ksi)+T(etta+1,gamma,ksi))/delX^2 +
...
(T(etta,gamma-1,ksi)-2*T(etta,gamma,ksi)+T(etta,gamma+1,ksi))/delY^2
+
(T(etta,gamma,ksi-1)-2*T(etta,gamma,ksi)+T(etta,gamma,ksi+1))/delZ2^2);etta
= 2:Nx-1gamma = 2:Ny-1ksi = 10:14_etta_gamma_ksi(etta,gamma,ksi) =
a3*delt*((T(etta-1,gamma,ksi)-2*T(etta,gamma,ksi)+T(etta+1,gamma,ksi))/delX^2 +
...
(T(etta,gamma-1,ksi)-2*T(etta,gamma,ksi)+T(etta,gamma+1,ksi))/delY^2
+
(T(etta,gamma,ksi-1)-2*T(etta,gamma,ksi)+T(etta,gamma,ksi+1))/delZ3^2);
end
Граничные условия при переходе из среды в среду
for etta = 2:Nx-1gamma = 2:Ny-1(etta,gamma,7) =
(lambtp1*T(etta,gamma,6)+lambtp2*T(etta,gamma,8))/(lambtp1+lambtp2);(etta,gamma,9)
= (lambtp1*T(etta,gamma,10)+lambtp2*T(etta,gamma,8))/(lambtp1+lambtp2);ksi =
2:Nz-1gamma = 2:Ny-1etta = 2:Nx-1(etta,gamma,ksi) =
T(etta,gamma,ksi)+delT_etta_gamma_ksi(etta,gamma,ksi);
Входные воздействияetta = 2:Nx-1gamma = 2:Ny-1
% T(etta,gamma,1) = delZ*utp/lambtp1 +
T(etta,gamma,2);(etta,gamma,1) = delZ*U(etta,gamma)/lambtp1 + T(etta,gamma,2);
Граничные условия:etta = 1:Nxgamma =1:Nyksi =
1:Nz(Nx,gamma,ksi) = T(Nx-1,gamma,ksi);(etta,Ny,ksi) =
T(etta,Ny-1,ksi);(etta,gamma,Nz) = T(etta,gamma,Nz-1);(Nx,gamma) =
T(Nx-1,gamma);(etta,Ny) = T(etta,Ny-1);etta = 1:Nxgamma = 1:Ny(etta,gamma) =
F(etta,gamma);(lmb,:) = T(4,2,6);(1), plot(lmb,T(4,4,12)), hold on(2),
plot(lmb,U), hold on
Расчет характеристик:
n = length(ForPlotT);
M = ForPlotT(n);_vv = U(4,4);= M/M_vv_x = pi*(i-0.5)/Lx;_y =
pi*(i-0.5)/Ly;= Psi_x^2 + Psi_y^2