УП_Лабы_Оптимизация управления ТП
.pdfПриложение В
Определение коэффициентов дифференциального уравнения по экспериментальной кривой разгона для объектов управления без самовыравнивания (часть 1, лабораторная работа № 6)
Рассмотрим приближенное решение дифференциального уравнения 30,433Y ′′(t)+ 3,822Y ′(t)= X (t) методом Рунге-Кутты. Для совпадения с траекторией, полученной по методу Эйлера, необходимо входное воздействие х (увеличение расхода металла из промежуточного ковша в кристаллизатор) задавать в кг/с, в нашем случае х=4.17кг/с. Количество шагов расчета n может совпадать с количеством шагов по методу Эйлера. Текст программы написан на языке PASCAL студентами гр. АМ-08:
program_k_dif_ur; uses crt,graph;
(*Инициализация массивов и переменных*) type mas=array[0..1000] of real; masy=array[0..1000] of real;
var h,a0,a1,a2,k1,k2,k3,k4,dy,dz,x :real; z,z0,y,y0 :real;
t,t1,ht,hy,max_y,min_y :real; sum1,sum2,integ1,integ2 :real; mas_t,mas_y :mas; y_o,t_o,dt,dy_o,k5,k6,k7 :masy; i,n,gd,gm,l :integer;
(*Инициализация функции расчета значения f*) function f(y,z:real):real;
begin f:=(x-a1*z-а0*y)/a2; end;
(*Ввод исходных данных для расчета*) begin
clrscr;
write('Введите количество шагов n :'); readln(n);
write('Введите коэффициент а1 :'); readln(а1);
write('Введите коэффициент а2 :'); readln(а2);
write('Введите входное воздействие x:');
172
readln(x);
write('Введите шаг по времени h :'); readln(h);
(*Расчет значений параметра выходной величины *) begin
t1:=0;
t:=t1;
y0:=0;
y:=y0;
a0:=0;
mas_y[0]:=0; mas_t[0]:=0; z0:=0; z:=z0;
l:=0;
for i:=1 to n do begin
if (l=20) then begin l:=0;
readln;clrscr;
end;
l:=l+1;
t:=t+h; y:=y0; z:=z0; k1:=h*f(y,z);
y:=y0+z0*h/2+k1*h/8; z:=z0+k1/2; k2:=h*f(y,z);
z:=z0+k2/2;
k3:=h*f(y,z); y:=y0+h*z0+h*k3/2; z:=z0+k3; dy:=h*(z0+(k1+k2+k3)/6); dz:=(k1+2*k2+2*k3+h*f(y,z))/6; y:=y0+dy; z:=z0+dz;
writeln('t=',t:5:3,' y(t)=', y:5:6); y0:=y; z0:=z; mas_y[i]:=y;mas_t[i]:=t;
end;
readln;
end;end.
Графики расчетных траекторий представлены на рис. В.1. Траектория, рассчитанная методом Эйлера (см. таблицу 1.3), практически совпадает с траекторией Z*(t) и на рисунке не отражена.
173
Рис. В.1. Расчетные траектории изменения уровня металла: Y(t) – действительное изменение уровня, Z(t) – по регистрирующему прибору методом Рунге-Кутта, Z*(t) – по регистрирующему прибору
точное решение уравнения
174
Приложение Г
Моделирование переходного процесса в системе с двухпозиционным регулированием (часть 2, лабораторная работа № 1)
Текст программы написан на языке VBA – Visual Basic for Applications.
Sub Двухпозиционный_регулятор()
Dim x As Double, y As Double
Dim z1(0 To 5000), z2(0 To 5000)
Dim dz1 As Double, dz2 As Double
Dim t As Double
Dim Tob As Integer |
'постоянная времени объекта |
Const tz = 10 |
'время запаздывания объекта |
Const Zzad = 180 |
'заданное значение температуры |
Const zn = 20 |
'зона нечувствительности регулятора |
y = 330 Tob = 450 z1(0) = 0 z2(0) = 0 dz1 = 0 dz2 = 0
Worksheets(1).Activate Cells(1, 1).Value = "t" Cells(1, 2).Value = "y"
Cells(1, 3).Value = "dz1" Cells(1, 4).Value = "z1" Cells(1, 5).Value = "dz2" Cells(1, 6).Value = "z2" Cells(2, 1).Value = 0 Cells(2, 2).Value = y Cells(2, 3).Value = dz1 Cells(2, 4).Value = z1(0) Cells(2, 5).Value = dz2 Cells(2, 6).Value = z2(0)
For t = 1 To 2000
z1(t) = dz1 + z1(t - 1)
dz1 = (1 / Tob) * (y - z1(t)) z2(t) = dz2 + z2(t - 1)
175
dz2 = (1 / tz) * (z1(t) - z2(t)) |
|
If (z2(t) <= (Zzad - zn / 2)) And (y = 0) Then |
|
y = 330 |
'нагрев |
Tob = 450 |
|
ElseIf (z2(t) > (Zzad + zn / 2)) And (y = 330) Then |
|
y = 0 |
'охлаждение |
Tob = 480
End If
Cells(t + 2, 1).Value = t
Cells(t + 2, 2).Value = y
Cells(t + 2, 3).Value = dz1
Cells(t + 2, 4).Value = z1(t)
Cells(t + 2, 5).Value = dz2
Cells(t + 2, 6).Value = z2(t)
Next t
End Sub
Температура,град
200
180
160
140
120
100
80
60
40
20
0
0 |
500 |
1000 |
1500 |
Время,сек
Рис. Г.1. Переходный процесс в системе двухпозиционного регулирования температуры
176
Приложение Д
Моделирование переходного процесса в системе с трехпозиционным регулированием (часть 2, лабораторная работа № 2)
Текст программы написан на языке VBA – Visual Basic for Applications.
Const a2 |
= 0.027 |
'коэффициенты |
Const a1 |
= 4.02 |
'статической |
Const a0 |
= 356.7 |
'характеристики |
'Линия регрессии
Function Yx(x As Double) As Double Yx = a0 + a1 * x + a2 * x ^ 2
End Function
Sub трехпозиционный_регулятор()
Dim x As Double, y As Double
Dim z1(0 To 5000), z2(0 To 5000)
Dim dz1 As Double, dz2 As Double
Dim sigma As Integer |
|
Dim t As Double |
|
Const Tob = 100 |
'постоянная времени объекта |
Const tz = 10 |
'время запаздывания объекта |
Const ki = 0.2 |
'скорость ИМ |
Const Zzad = 670 |
'заданное значение температуры |
Const zn = 30 |
'зона нечувствительности регулятора |
x = 46 |
'начальное положение вала ИМ |
sigma = 1 |
|
z1(0) = Yx(x) |
|
z2(0) = Yx(x) |
|
dz1 = 0 |
|
dz2 = 0 |
|
Worksheets(1).Activate |
|
Cells(1, 1).Value = "t" |
|
Cells(1, 2).Value = "x" |
|
Cells(1, 3).Value = "y" |
|
Cells(1, 4).Value = "dz1" |
|
Cells(1, 5).Value = "z1" |
|
Cells(1, 6).Value = "dz2" |
|
Cells(1, 7).Value = "z2" |
|
Cells(1, 8).Value = "sigma" |
|
177
Cells(2, 1).Value = 0
Cells(2, 2).Value = x
Cells(2, 3).Value = Yx(x)
Cells(2, 4).Value = dz1
Cells(2, 5).Value = z1(0)
Cells(2, 6).Value = dz2
Cells(2, 7).Value = z2(0)
Cells(2, 8).Value = sigma
For t = 1 To 2000 |
|
If x = 0 Then x = x |
'ограничение на |
If x = 100 Then x = x |
'концевые выключатели ИМ |
x = x + ki * sigma |
|
y = Yx(x) |
|
z1(t) = dz1 + z1(t - 1) |
|
dz1 = (1 / Tob) * (y - z1(t)) |
|
z2(t) = dz2 + z2(t - 1) |
|
dz2 = (1 / tz) * (z1(t) - z2(t)) |
|
Cells(t + 2, 1).Value = t |
|
Cells(t + 2, 2).Value = x |
|
Cells(t + 2, 3).Value = y |
|
Cells(t + 2, 4).Value = dz1 |
|
Cells(t + 2, 5).Value = z1(t) |
|
Cells(t + 2, 6).Value = dz2 |
|
Cells(t + 2, 7).Value = z2(t) |
|
Cells(t + 2, 8).Value = sigma |
|
If Zzad - z2(t) > zn / 2 Then sigma = 1
If Zzad - z2(t) > -zn / 2 And Zzad - z2(t) < zn / 2 Then sigma = 0 If Zzad - z2(t) < -zn / 2 Then sigma = -1
Next t
End Sub
178
|
720 |
|
|
75 |
|
|
700 |
|
|
70 |
|
Температура,град |
680 |
|
|
65 |
Положение вала ИМ,% |
660 |
|
|
60 |
||
640 |
|
|
55 |
||
620 |
|
|
50 |
||
|
|
|
|
||
|
600 |
|
|
45 |
|
|
580 |
|
|
40 |
|
|
0 |
500 |
1000 |
1500 |
|
|
|
|
Время,сек |
|
|
Рис. Д.1. Переходный процесс в системе трехпозиционного |
|||||
|
|
регулирования температуры |
|
|
|
179
Приложение Е
Моделирование переходного процесса в системе с ПИ-регулятором (часть 2, лабораторная работа № 3)
Текст программы написан на языке VBA – Visual Basic for Applications.
Const a2 = -0.0269 'коэффициенты
Const a1 = 9.88 'статической
Const a0 = 49.16 'характеристики
'Линия регрессии
Function Yx(x As Double) As Double Yx = a0 + a1 * x + a2 * x ^ 2
End Function
Sub ПИ_регулятор()
Dim Xd As Double, Xp As Double, y As Double
Dim z1(0 To 5000), z2(0 To 5000)
Dim dz1 As Double, dz2 As Double
Dim sigma As Integer
Dim epsilon As Double, integral As Double
Dim t As Double, i As Double |
|
Const Tob = 75 |
'постоянная времени объекта |
Const tz = 5 |
'время запаздывания объекта |
Const Tiz = 50 |
|
Const kp = 0.7 |
|
Const ki = 0.83 |
'скорость ИМ |
Const Zzad = 550 |
'заданное значение температуры |
Const zn = 2 |
'зона нечувствительности регулятора |
Xd = 47 |
'начальное положение вала ИМ |
Xp = 0 |
|
sigma = 1 |
|
integral = 0 |
|
z1(0) = Yx(Xd) |
|
z2(0) = Yx(Xd) |
|
dz1 = 0 |
|
dz2 = 0 |
|
Worksheets(1).Activate |
|
Cells(1, 1).Value = "t" |
|
180
Cells(1, 2).Value = "Xd"
Cells(1, 3).Value = "y"
Cells(1, 4).Value = "dz1"
Cells(1, 5).Value = "z1"
Cells(1, 6).Value = "dz2"
Cells(1, 7).Value = "z2"
Cells(1, 8).Value = "sigma"
Cells(1, 9).Value = "epsilon"
Cells(1, 10).Value = "integral"
Cells(1, 11).Value = "Xp"
Cells(2, 1).Value = 0
Cells(2, 2).Value = Xd
Cells(2, 3).Value = Yx(Xd)
Cells(2, 4).Value = dz1
Cells(2, 5).Value = z1(0)
Cells(2, 6).Value = dz2
Cells(2, 7).Value = z2(0)
Cells(2, 8).Value = sigma
Cells(2, 9).Value = epsilon
Cells(2, 10).Value = integral
Cells(2, 11).Value = Xp
For t = 1 To 500
epsilon = Zzad - z2(t - 1) If Xp >= 100 Then integral = integral
Xp = 100 Else
integral = integral + epsilon
Xp = kp * (epsilon + (1 / Tiz) * integral) End If
Xd = Xd + ki * sigma y = Yx(Xd)
z1(t) = dz1 + z1(t - 1)
dz1 = (1 / Tob) * (y - z1(t)) z2(t) = dz2 + z2(t - 1)
dz2 = (1 / tz) * (z1(t) - z2(t)) If Xp - Xd > 0 Then
If Xp - Xd + zn > 0 Then sigma = 1
End If
If Xp - Xd + zn < 0 Then sigma = 0
181