Определение максимального шага интегрирования

Для определения максимального шага интегрирования введем критерий характеризующий качество сохранения процесса. В качестве данного критерия используем величину математического ожидания ошибки воспроизведения .

Для определения значения критерия воспользуемся экспортом данных с самописца (scope в среде Simulink) в рабочую среду MatLab, в перемену с именем s_ur в виде массива значений (array). И командами:

sim('m');

sim('m')

i_ur = 0;

for i = 1:size(s_ur)-1

i_ur(i) = abs(s_ur(i,2)-s_ur_i(round(s_ur(i,1)*100000+1), 2));

i = i + 1;

end

В результате выполнения i_ur содержит массив ошибки воспроизведения . Для вычисления математического ожидания воспользуемся командой mean()

В качестве значения максимального шага интегрирования выбирается максимальное значении шага при котором значение mean(i_ur) не превышает одну сотую. Для поиска максимального шага будем использовать только две значащие цифры шага. Результаты поиска представлены в таблице 2.

таблица 2

Название метода	Шаг	mean(i_ur)
Эйлера	0.0012	0.8864
0.0013	0.9706
0.0014	1.0432
Рунге-Кутта 2-го порядка	0.015	0.6919
0.016	0.9063
0.017	1.0199
Рунге-Кутта 3-го порядка	0.020	0.4708
0.021	0.8764
0.022	1.0064
Рунге-Кутта 4-го порядка	0.024	0.4063
0.025	0.9064
0.026	1.1062
Рунге-Кутта 5-го порядка	0.026	0.4190
0.027	0.8765
0.028	1.0114
Адамса	0.009	0.7668
0.010	0.9548
0.011	1.4445
Гира	0.0021	0.8767
0.0022	0.9992
0.0023	1.1780

В таблице, для многошаговых методов Адамса и Гира, под шагом подразумевается максимальный шаг (Max step size).

Выбор метода для дальнейшего интегрирования необходимо проводить с учетом реального времени, затрачиваемого на моделирование.

Для измерения этого времени воспользуемся следующей конструкцией:

tic; sim(‘m’); toc

Время моделирования 5000 секунд процесса различными методами с использованием максимального шага интегрирования приведено в таблице 3.

таблица 3

Название метода	Шаг	время моделирования в секундах
Эйлера	0.0013
Рунге-Кутта 2-го порядка	0.016	3.9
Рунге-Кутта 3-го порядка	0.021	3.2
Рунге-Кутта 4-го порядка	0.025	3.6
Рунге-Кутта 5-го порядка	0.027	5.6
Адамса	0.010
Гира	0.0022

В качестве основного метода интегрирования для дальнейшего исследования был выбран метод Рунге-Кутта 3-го порядка с шагом 0.021.

Выбор данного метода в качестве основного обоснован минимальным временем моделирования пяти тысяч секунд процесса при сохранении качества процессов.

 

Синтез «в большом»

За показатель качества примем:

,

где -максимальное относительное отклонение от установившегося значения,

время вхождения в 1% зону

 

Для изучения примем первые 35 секунд процесса (скачкообразное изменение с 40.79 на 35 происходит на первой секунде). В качестве наблюдаемой величины воспользуемся .

Для определения начального вектора приближения:

В качестве примем минимальное значение при котором, при процесс автоколебательный.

,

рис. 14

,

рис. 15

Для определения значения параметра при заданном перейдем к рассмотрению более мелкого участка, а именно первых 10 секунд процесса.

Найдём минимальное значение , при котором не превышает 7 секунд.

,

рис. 16

Примем максимально допустимое отклонение в 10%, тогда

,

 

рис. 17

 

Начальный вектор имеет следующее значение: ,

Первая итерация поиска представлена в таблицах 4 и 5.

таблица 4

		0.1	0.15	0.2	0.25	0.3
2.15	14,3131%-9,97c	10,1862%-10c	9,4692%-6,462c	8,8826%-5,03c	8,3924%-3,998c
3.225	11,5734%-9,524c	10,5275%-6,076c	9,719%-3,811c	9,0723%-3,488c	8,5408%-3,263c
4.3	11,8732%-7,141c	10,7278%-4,223c	9,8608%-3,726c	9,1774%-3,445c	8,6216%-3,273c
5.375	12,0786%-4,926c	10,8605%-4,079c	9,9527%-3,663c	9,2445%-3,479c	8,6725%-3,333c
6.45	12,2292%-4,851c	10,9552%-3,146c	10,0172%-3,706c	9,291%-3,558c	8,7076%-3,414c

таблица 5

		0.1	0.15	0.2	0.25	0.3
2.15	11,27	10,06	7,36	6,19	5,32
3.225	10,14	7,41	5,58	5,16	4,85
4.3	8,56	6,17	5,57	5,16	4,88
5.375	7,07	6,11	5,55	5,21	4,93
6.45	7,06	5,49	5,60	5,28	5,00

В качестве нового вектора для поиска оптимального примем: ,

Вторая итерация поиска представлена в таблицах 6 и 7.

таблица 6

		0,15	0,23	0,30	0,38	0,45
1,6	14,3131%-9,97c	10,1862%-10c	9,4692%-6,462c	8,8826%-5,03c	8,3924%-3,998c
2,4	11,5734%-9,524c	10,5275%-6,076c	9,719%-3,811c	9,0723%-3,488c	8,5408%-3,263c
3,25	11,8732%-7,141c	10,7278%-4,223c	9,8608%-3,726c	9,1774%-3,445c	8,6216%-3,273c
4,06	12,0786%-4,926c	10,8605%-4,079c	9,9527%-3,663c	9,2445%-3,479c	8,6725%-3,333c
4,9	12,2292%-4,851c	10,9552%-3,146c	10,0172%-3,706c	9,291%-3,558c	8,7076%-3,414c

таблица 7

		0,15	0,23	0,30	0,38	0,45
1,6	10,70	8,72	6,29	4,89	4,58
2,4	9,51	5,87	5,23	4,46	4,21
3,25	6,83	5,32	4,85	4,48	4,24
4,06	6,17	5,31	4,87	4,51	4,27
4,9	6,16	5,31	4,91	4,89	4,30

В качестве нового вектора для поиска оптимального примем

,

Необходимо учесть, что не может принимать значения выходящие за пределы [-10;10] В во время работы. Уже при значение при переходе с номинального на заданный режим изменяется в диапазоне от 5.3 до 8 В. Дальнейшее увеличение считаю не допустимым, а в качестве оптимального принимаю текущий вектор параметров. А именно:

,

Графики переходных процессов при данной настройке регулятора приведены на рисунке 18

рис. 18

 

Линеаризация

Для найденного равновесного состояния произведем линеаризацию системы. После упрощения имеем следующую передаточную функцию:

Структурная схема линеаризованной системы представлена на рисунке 19.

рис. 19

 

Синтез «в малом»

По линеаризованной модели системы найдем параметры регулятора, обеспечивающие приемлемое качество процессов в окрестности равновесного состояния. Для поиска оптимального регулятора воспользуемся методом аналогичным примененному в пункте синтез «в большом».

Показателя качества:

, где – перерегулирование и - время регулирования

В качестве вектора начального приближения примем полученный ранее, а именно:

,

Первая итерация поиска представлена в таблицах 8 и 9.

таблица 8

		0,23	0,34	0,45	0,56	0,68
1,20	84,259%-73,2616 c	70,3943%-10,2246c	62,3809%-6,1692c	59,1957%-4,8913c	56,1342%-3,8685c
1,80	62,863%-10,3183c	53,3107%-5,3883c	49,0721%-3,9875c	47,6026%-3,5419c	47,5395%-3,2098c
2,40	50,459%-5,6059c	43,9957%-4,4272c	42,1346%-3,2639c	42,503%-3,4272c	43,1786%-2,7044c
3,00	41,9152%-4,5619c	36,1826%-3,8033c	36,8609%-3,2613c	39,06%-2,9174c	40,2409%-2,7222c
3,60	34,4102%-4,4723c	32,2137%-3,6853c	33,4198%-3,2837c	36,2158%-2,9601c	38,5755%-2,7344c

таблица 9

		0,23	0,34	0,45	0,56	0,68
1,20	76,56	28,28	23,03	21,18	19,55
1,80	26,08	19,77	17,51	16,76	16,51
2,40	19,06	16,30	14,93	15,15	14,85
3,00	15,77	13,52	13,34	13,76	13,98
3,60	13,45	12,24	12,32	12,94	13,49

В качестве нового вектора для поиска оптимального примем: ,

 

Вторая итерация поиска представлена в таблицах 10 и 11.

таблица 10

		0,17	0,26	0,34	0,43	0,51
1,80	73,6711%-21,0125c	59,3017%-8,0054c	53,3107%-5,3883c	49,5809%-4,0899c	48,0388%-3,722c
2,70	53,0305%-7,5487c	43,7773%-4,3253c	40,331%-3,7183c	39,4406%-3,3056c	39,9355%-3,0534c
3,60	40,1934%-5,3366c	33,2308%-4,1432c	32,2137%-3,6853c	33,083%-3,3478c	34,7238%-3,1096c
4,50	31,8629%-4,245c	28,0157%-3,2891c	28,2507%-2,7891c	30,466%-3,3061c	31,8963%-3,1714c
5,40	25,5611%-3,9338c	25,088%-3,305c	25,7281%-2,2752c	27,8785%-3,4433c	30,496%-3,2073c

таблица 11

		0,17	0,26	0,34	0,43	0,51
1,80	36,81	23,39	19,77	17,74	17,02
2,70	21,19	16,16	14,70	14,15	14,12
3,60	15,79	12,87	12,24	12,27	12,59
4,50	12,53	10,71	10,43	11,45	11,79
5,40	10,42	9,84	9,31	10,77	11,39

В качестве нового вектора для поиска оптимального примем: ,

Третья, четвертая, пятая, шестая итерации выполняются аналогичным образом и здесь не приводятся.

Седьмая итерация поиска представлена в таблицах 12 и 13.

таблица 12

		0,03	0,04	0,05	0,06	0,08
9,50	35,5616%-10,7122c	31,018%-8,0145c	26,1047%-7,1982c	22,2173%-5,3879c	16,1047%-4,6272c
14,24	19,0072%-7,0915c	13,6959%-6,2795c	9,9657%-4,1794c	7,3857%-3,8057c	2,5153%-2,4512c
18,99	9,583%-5,626c	4,4244%-4,7035c	0,5161%-2,5545c	0%-2,9278c	0%-5,5828c
23,74	2,1337%-4,7655c	0%-3,603c	0%-5,2288c	0%-6,0684c	0%-7,4232c
28,49	0%-4,5358c	0%-6,3626c	0%-6,9665c	0%-7,8558c	0%-9,3798c

таблица 13

		0,03	0,04	0,05	0,06	0,08
9,50	18,17	14,92	12,87	10,44	8,07
14,24	10,67	8,50	5,92	4,88	2,47
18,99	6,81	4,62	1,94	2,05	3,91
23,74	3,98	2,52	3,66	4,25	5,20
28,49	3,18	4,45	4,88	5,50	6,57

Т.к. изменение значения параметров регулятора не увеличивает качество процесса (показатель качества в «центре»), то в качестве вектора параметров оптимального регулятора примем: , .

Т.о. передаточная функция найденного регулятора имеет вид:

Переходный процесс представлен на рисунке 20.

рис. 20

Время регулирования 2,5545c

Перерегулирование 0,5161%