Помощь в написании студенческих работ
Антистрессовый сервис

Расчет ПИД и ПИ-регулятора

КурсоваяПомощь в написанииУзнать стоимостьмоей работы

Так как выходной сигнал имеет конечное установившееся значение, то есть система приходит к статическому режиму, в котором скорости изменения входного и выходного сигналов равны нулю, то можно говорить о том, что объект с самовыравниванием. Для построения амплитудно-фазовой характеристики рабочей модели достаточно в её передаточной функции сделать подстановку S = j· щ и теперь, при изменении… Читать ещё >

Расчет ПИД и ПИ-регулятора (реферат, курсовая, диплом, контрольная)

КУРСОВАЯ РАБОТА

Расчет ПИД и ПИ регулятора

Цель работы Задание Исходные данные

1. Построение переходной кривой объекта

2. По переходной кривой методом «площадей» Симою М. П. определить параметры нескольких моделей объекта (площадь S рассчитать вручную)

3. Расчёт переходной кривой по передаточной функции. Выбор рабочей модели.

4. Построение АФХ рабочей модели

5. Выбор закона регулирования

6. Построение области устойчивости в плоскости настроечных параметров регулятора

7. Расчет и построение в плоскости параметров настроек кривой равного значения

8. Определение оптимальных параметров регулятора

9. Построение АФХ разомкнутой и АЧХ замкнутой АСР по задающему воздействию для оптимальных настроек регулятора

9.1. Передаточные функции разомкнутых систем

9.1.1 ПИ-регулятор

9.1.2 ПИД-регулятор

9.2 Построение АФХ

9.2.1 ПИ-регулятор

9.2.2 ПД-регулятор

10. Построение переходных кривых в замкнутой АСР по задающему и возмущающему воздействию методом Акульшина

11. Анализ качества регулирования

11.1 Статическая ошибка еуст

11.2 Время регулирования

11.2.1 Для кривой по задающему воздействию

11.2.2 Для кривой по возмущающему воздействию

11.3 Перерегулирование

11.4 Степень затухания

11.4.1 Для кривой по задающему воздействию

11.4.2 Для кривой по возмущающему воздействию

12. Выбор промышленного регулятора

13. Список используемой литературы

Цель работы

Определить настройки типового регулятора (ПИ, ПИД), минимизирующие интегральный квадратичный критерий I0 при заданном ограничении Mзад (Mе зад) или m>mзад. Выбрать промышленный регулятор и его настройки.

Задание

1. Построить переходную кривую объекта по табличным данным.

2. По переходной кривой методом «площадей» Симою М. П. определить параметры нескольких моделей объекта (площадь S1 рассчитать вручную).

3. По найденным передаточным функциям методом обратного преобразования Лапласа рассчитать и построить переходные кривые моделей (две точки одной из кривых рассчитать вручную). Выбрать рабочую модель, наиболее близкую к объекту.

4. Построить нормальную и расширенную АФХ рабочей модели объекта (одну точку АФХ рассчитать вручную).

5. Выбрать закон регулирования (расчет вести для двух законов регулирования).

6. Построить область устойчивости в плоскости настроечных параметров регулятора (одну точку кривой Д-разбиения для одного из регуляторов построить вручную).

7. Рассчитать и построить в плоскости параметров настроек кривую равного значения показателя колебательности М = Мзад = 1.5.

8. Определить оптимальные параметры регулятора.

9. Построить АФХ разомкнутой АСР (одну точку рассчитать вручную) и АЧХ замкнутой по задающему воздействию для оптимальных настроек регулятора.

10. Построить переходные кривые в замкнутой АСР по задающему и возмущающему воздействию методом Акульшина.

11. Провести анализ качества регулирования. Выбрать наилучший закон регулирования.

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

расчет параметр регулятор пи пид

Исходные данные

x = 20 кПа — амплитуда входного сигнала;

yуст = 30 °C — выходная (регулируемая) величина;

зап = 1 мин — запаздывание;

Tшк = 600 °C — диапазон шкалы;

Мзад = 1.5 — степень колебательности.

Таблица 1. Переходный процесс объекта:

t, мин

Дy, 0C

1. Построение переходной кривой объекта

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

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

По данным (Таблица 1) строится переходная кривая объекта (Рисунок 1), при этом запаздывание не учитывается.

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

Рисунок 1 — Переходная кривая объекта

2. Определение параметров нескольких моделей объекта по переходной кривой методом «площадей» Симою М. П

Математической моделью называется система математических соотношений (уравнений), устанавливающих связь между входными и выходными сигналами объекта.

В данном случае общий вид модели будет следующий:

— нормированная передаточная функция;

— коэффициент усиления ;

— время запаздывания (по исходным данным мин);

Нормированной передаточной функции соответствует нормированная переходная характеристика (t), которая определяется как отношение текущего значения выходного сигнала к его установившемуся значению: .

Для определения коэффициентов ai и bi нормированной передаточной функции используется метод «площадей» Симою.

(*)

sk - «площади» Симою; вычисляются по переходной кривой.

При известных «площадях» Симою, задаваясь определённой структурой модели можно определить её параметры (коэффициенты). «Площади» Симою определяются с помощью вспомогательной (t) функции:

.

(**)

; ;

— моменты вспомогательной функции.

Если из выражения (**) выразить, а затем приравнять правые части уравнений (*) и (**), то легко найти связь между моментами вспомогательной функции и «площадями» Симою:

Так — площадь под кривой вспомогательной функции.

Для расчёта площади s1 необходимо рассчитать значения вспомогательной функции (Таблица 2).

Таблица 2. Результаты расчёта вспомогательной функции

t, min

Y (t)

()

()

0,00

1,00

0,07

0,93

0,17

0,83

0,27

0,73

0,33

0,67

0,40

0,60

0,50

0,50

0,60

0,40

0,83

0,17

0,90

0,10

1,00

0,00

1,00

0,00

1,00

0,00

По данным Таблицы 2 строится график вспомогательной функции (Рисунок 2). Методом трапеций определяется площадь под кривой вспомогательной функции:

где t = 1 мин — шаг по времени.

Полученное значение и есть значение «площади» Симою s1.

Остальные расчёты проведём на ЭВМ (программа Simou. exe)

Рисунок 2 — Вспомогательная функция f (t)

РЕЗУЛЬТАТЫ РАСЧЕТА:

коэффициент усиления переходной функции KY = 1.5

Моменты

M0 = 5.466

M1 = -18.888

M2 = 47.711

M3 = -94.925

M4 = 156.064

Площади Симою:

S0 = 1

S1 = 5.466

S2 = 10.995

S3 = 4.560

S4 = -16.866

S5 = -16.605

Получил параметры 5 моделей.

Параметры передаточной функции модели № 1:

Коэффициент усиления K = 1,5

Запаздывание: фзап = 0

Коэффициент числителя (степень m = 0): b[0] = 1

Коэффициент знаменателя (степень n = 1): a[0] = 1 a[1] = 5,46

Вид передаточной функции:

Параметры передаточной функции модели № 2:

Коэффициент усиления: K = 30

Запаздывание: фзап = 0

Коэффициент числителя (степень m = 0): b[0] = 1

Коэффициент знаменателя (степень n = 1): a[0] = 5.466, a[1] = 10.995

Вид передаточной функции:

Параметры передаточной функции модели № 3:

Коэффициент усиления: K = 30

Запаздывание: фзап = 0

Коэффициент числителя (степень m = 0): b[0] = 1

Коэффициент знаменателя (степень n = 3): a[1] = 5.466, a[2] = 10.995, a[3] = 4.560

Вид передаточной функции:

Параметры передаточной функции модели № 6:

Коэффициент усиления: K = 30

Запаздывание: фзап = 0

Коэффициент числителя (степень m = 1): b[0] = 1, b[1] = -0.414

Коэффициент знаменателя (степень n = 2): a[0] = 1, a[1] = 5.051, a[2] = 8.728

Вид передаточной функции:

Параметры передаточной функции модели № 7:

Коэффициент усиления: K = 30

Запаздывание: фзап = 0

Коэффициент числителя (степень m = 1): b[0] = 1, b[1] = 3.697

Коэффициент знаменателя (степень n = 3): a[0] = 1, a[1] = 9.164, a[2] = 31.211, a[3] = 45.222

Вид передаточной функции:

Расчет переходных кривых проведем для вариантов 2, 3, 7, т.к. только для них выполняется критерий устойчивости Стодолы.

3. Расчёт переходной кривой по передаточной функции. Выбор рабочей модели

Расчёт переходной кривой по передаточной функции для варианта 3:

Корни A (s):

Как видно из рис. 3 переходная кривая третий вспомогательной модели объекта наиболее точно совпадает с переходной кривой самого объекта.

Передаточная функция данной модели:

.

Рисунок 3 — Сравнение переходных кривых

4. Расчет АФХ рабочей модели объекта. Построение нормальной АФХ модели

Для построения амплитудно-фазовой характеристики рабочей модели достаточно в её передаточной функции сделать подстановку S = j· щ и теперь, при изменении частоты от 0 до ?, комплексный вектор WM(j· щ) будет поворачиваться, описывая годограф АФХ.

В

ручную рассчитаем одну точку АФХ:

Результаты полностью совпадают с данными, полученными на ЭВМ, что говорит о правильности найденного выражения для построения АФХ.

ПАРАМЕТРЫ ПЕРЕДАТОЧНОЙ ФУНКЦИИ

Степень числителя m = 0

Полином числителя: b0 = 1.00

Степень знаменателя n = 3

Полином знаменателя: a0 = 1.00

a1 = 5.46

a2 = 10.99

a3 = 4.56

Коэффициент усиления Ку = 1.50

Запаздывание tau = 1.00

ПАРАМЕТРЫ РАСЧЕТА

Начальная частота fo 0.00

Число точек n 20

Шаг по частоте h 0.10

РЕЗУЛЬТАТЫ РАСЧЕТА Таблица 3. Результаты расчёта нормальной АФХ рабочей модели Рисунок 4 — Нормальная АФХ рабочей модели объекта

5. Выбор закона регулирования

Пропорционально-интегральный (ПИ) регулятор (WПИ = k1 + k0 / S)

Уравнение ПИ-регулятора во временной области:

Передаточная и переходная функции ПИ-регулятора:

Пропорционально-дифференциальный (ПД) регулятор (WПД = kр + Tд?S)

Уравнение ПД-регулятора во временной области:

Регулятор является параллельным соединением П и Д регуляторов имеет два параметра настройки, и .

Передаточная и переходная характеристики регулятора

.

при t?0, где — дельта-функция,

6. Построение области устойчивости в плоскости настроечных параметров регулятора

Кривая D-разбиения является границей области устойчивости и показывает область изменения настроечных параметров регулятора, при которых система устойчива.

Кривая D-разбиения может быть получена из характеристического уравнения замкнутой АСР подстановкой s = j:

W (s) + 1 = 0, что эквивалентно Dз(s) = 0;

Передаточная функция разомкнутой АСР:

W (s) = Wр(s)· Wм(s), где Wр(s) — передаточная функция регулятора.

Уравнение границы области устойчивости:

Wр(s)· Wм(s) + 1 = 0.

ПИ-регулятор:

Преобразуем это уравнение следующим образом

K1•V (s) + K0.X (s) + 1 = 0, где V (s) = Wм(s); X (s) = Wм(s)/s;

Сделаем подстановку s = j: K1· V (j) + K0· X (j) + 1 = 0.

Получим систему уравнений:

K1V1() + K0X1() +1 = 0, где V1() = Re V (j); X1()= Re X (j).

K1V2() + K0X2() = 0.

V2() = Im V (j); X2() = Im X (j)

Решение системы можно найти методом определителей:

К1 =, К0 = ,

.

Re () и Im () были получены в пункте 4.

Рассчитаем K1 и K0 для = 0.3:

Рассчитываем значение кривой D-разбиения в точке соответствующей частоте

V1() = Re Wм(j); V2()= Im Wм(j)

кПа/(0C*мин);

кПа/0C.

ИСХОДНЫЕ ДАННЫЕ ДЛЯ РАСЧЕТА Параметры передаточной функции:

Коэффициент усиления K =1,5

Запаздывание tau = 1

Коэффициенты числителя (степень m = 0):

b[0] = 1

Коэффициенты знаменателя (степень n = 3):

a[0] = 1 a[1] = 5,46 a[2] = 10,99 a[3] = 4,56

Расчет для АСР с ПИ-регулятором.

Расчет границы устойчивости.

РЕЗУЛЬТАТЫ РАСЧЕТА Кривая D-разбиения Таблица 5. Результат расчета границы устойчивости для ПИ-регулятора Рисунок 7 — Кривая D-разбиения границы устойчивости для ПИ-регулятора

ПД-регулятор:

ИСХОДНЫЕ ДАННЫЕ ДЛЯ РАСЧЕТА Параметры передаточной функции:

Коэффициент усиления K =1,5

Запаздывание tau = 1

Коэффициенты числителя (степень m = 0):

b[0] = 1

Коэффициенты знаменателя (степень n = 3):

a[0] = 1 a[1] = 5,46 a[2] = 10,99 a[3] = 4,56

Расчет границы устойчивости.

РЕЗУЛЬТАТЫ РАСЧЕТА Кривая D-разбиения в программе ТАУ 1.0 для DOS.

Таблица 6. Результат расчета границы устойчивости для ПД-регулятора.

Рисунок 8 — Кривая D-разбиения границы устойчивости для ПД-регулятора

7. Расчет и построение в плоскости параметров настроек кривую равного значения: показателя колебательности M = Mзад = 1.5

ПИ-регулятор:

ИСХОДНЫЕ ДАННЫЕ ДЛЯ РАСЧЕТА Параметры передаточной функции:

Коэффициент усиления K =1,5

Запаздывание tau = 1

Коэффициенты числителя (степень m = 0):

b[0] = 1

Коэффициенты знаменателя (степень n = 3):

a[0] = 1 a[1] = 5,46 a[2] = 10,99 a[3] = 4,56

Расчет для АСР с ПИ-регулятором.

Расчет кривой для заданного показателя колебательности Mзад = 1,5

РЕЗУЛЬТАТЫ РАСЧЕТА Кривая D-разбиения Таблица 7. Результат расчета границы устойчивости для ПИ-регулятора при заданном показателе колебательности Мзад = 1.5.

Рисунок 9 — Кривая D-разбиения для ПИ-регулятора при заданном показателе колебательности Мзад = 1.5.

ПД-регулятор:

ИСХОДНЫЕ ДАННЫЕ ДЛЯ РАСЧЕТА Параметры передаточной функции:

Коэффициент усиления K =1,5

Запаздывание tau = 1

Коэффициенты числителя (степень m = 0):

b[0] = 1

Коэффициенты знаменателя (степень n = 3):

a[0] = 1 a[1] = 5,46 a[2] = 10,99 a[3] = 4,56

Расчет для АСР с ПИ-регулятором.

Расчет кривой для заданного показателя колебательности Mзад = 1,5

РЕЗУЛЬТАТЫ РАСЧЕТА Кривая D-разбиения

Таблица 8. Результат расчета границы устойчивости для ПД-регулятора при заданном показателе колебательности Mзад = 1,5.

Рисунок 10 — Кривая D-разбиения для ПД-регулятора при заданном показателе колебательности Мзад = 1.5.

8. Определение оптимальных параметров регулятора

Оптимальные параметры будем выбирать из условия минимизации интегрального квадратичного критерия I0 на кривой М = Мзад:

Положение оптимальной (рабочей) точки, как в случае ПИ-, так и в случае ПД-регулятора существенно зависит от степени неопределенности задачи.

Известно возмущение и передаточная функция объекта по каналу возмущения. В нашем случае скачкообразное возмущение приложено со стороны регулирующего органа, а передаточная функция рассчитывается по переходной кривой. В этом случае рабочая точка лежит несколько правее максимума граничной кривой и определяется формулами:

щр = 1,2*щmax, [1/мин]

или щр = 0,67*щп, [1/мин]

Найдём оптимальные параметры регуляторов.

1) ПИ-регулятор:

Пользуясь таблицей и рисунком, можно определить:

щр = 1,2*щmax = 1,2*0,31 = 0,37 [1/мин]

К0 = 0,14 [кПа/(0С*мин)]

К1 = 0,82 [кПа/0С].

2) ПД-регулятор:

Пользуясь таблицей и рисунком, можно определить:

щр = 1,2*щmax = 1,2*0,75 = 0,9 [1/мин]

К1 = 1,61 [кПа/0С]

К2 = 3,75 [кПа/0С * мин].

9. Построение АФХ разомкнутой и АЧХ замкнутой АСР по задающему воздействию для оптимальных настроек регулятора

9.1 Передаточные функции разомкнутых систем

9.1.1 ПИ-регулятор

9.1.2 ПД-регулятор

9.2 Построение АФХ

Для расчёта АФХ в передаточную функцию разомкнутой системы сделаем подстановку s = j:

9.2.1 ПИ-регулятор

АФХ разомкнутой АСР с ПИ-регулятором:

ИСХОДНЫЕ ДАННЫЕ ДЛЯ РАСЧЕТА Параметры передаточной функции регулятора:

Коэффициент усиления K = 1

Запаздывание tau = 0

Коэффициенты числителя (степень m = 1):

b[0] = 0,14 b[1] = 0,82

Коэффициенты знаменателя (степень n = 1):

a[0] = 0 a[1] = 1

Параметры передаточной функции:

Коэффициент усиления K =1,5

Запаздывание tau = 1

Коэффициенты числителя (степень m = 0):

b[0] = 1

Коэффициенты знаменателя (степень n = 4):

a[0] = 1 a[1] = 5,46 a[2] = 10,99 a[3] = 4,56

РЕЗУЛЬТАТЫ РАСЧЕТА Обозначены: Re — действительная часть АФХ, Im — мнимая часть АФХ) Частотные характеристики разомкнутой системы Таблица 9. Результаты расчета АФХ разомкнутой системы с ПИ-регулятором.

Рисунок 13 — АФХ разомкнутой системы с ПИ-регулятором

АФХ разомкнутой АСР с ПД-регулятором:

ИСХОДНЫЕ ДАННЫЕ ДЛЯ РАСЧЕТА Параметры передаточной функции регулятора:

Коэффициент усиления K = 1

Запаздывание tau = 0

Коэффициенты числителя (степень m = 2): b[0] = 3,75 b[1] = 1.61

Коэффициенты знаменателя (степень n = 1): a[0] = 0

Параметры передаточной функции:

Коэффициент усиления K =1,5

Запаздывание tau = 1

Коэффициенты числителя (степень m = 0): b[0] = 1

Коэффициенты знаменателя (степень n = 3): a[0] = 1 a[1] = 5,46 a[2] = 10,99 a[3] = 4,56

РЕЗУЛЬТАТЫ РАСЧЕТА Обозначены: Re — действительная часть АФХ, Im — мнимая часть АФХ) Частотные характеристики разомкнутой системы Таблица 10. Результаты расчета АФХ разомкнутой системы с ПД-регулятором.

Рисунок 14 — АФХ разомкнутой системы с ПД-регулятором

АЧХ замкнутой АСР по задающему воздействию

Таблица 11. Частотные характеристики замкнутой системы по заданию ПИ-регулятор

ПД-регулятор

Рисунок 15 — АЧХ замкнутой АСР с ПИ-регулятором Рисунок 16 — АЧХ замкнутой АСР с ПД-регулятором

10. Построение переходных кривых в замкнутой АСР по задающему и возмущающему воздействию методом Акульшина

Анализ качества регулирования может производиться по прямым и косвенным критериям качества. Для анализа по прямым критериям необходимо построение переходного процесса. Методы построения переходных процессов делятся на две основные категории — точные и приближённые. Точные методы подразумевают точное решение дифференциальных уравнений, которыми описывается система. На практике чаще используются приближённые и в основном частотные методы.

Частотный метод Акульшина основан на том, что на вход АСР подаётся прямоугольная волна. Период выбирается так, что, то есть к моменту времени t = 0 процесс установится, и начальные условия можно считать нулевыми. Прямоугольная волна раскладывается в ряд, значит, на выходе системы также получаем сигнал, состоящий из нескольких гармоник (по принципу суперпозиции). Для расчёта достаточно 50 гармоник.

Частота прямоугольной волны принимается равной .

Таблица 12. Значение переходного процесса в АСР с ПИ-регулятором

T

НВ

0,00

0,00

0,00

3,00

0,14

0,15

6,00

0,76

0,69

9,00

1,18

0,79

12,00

1,13

0,41

15,00

0,92

0,05

18,00

0,85

— 0,03

21,00

0,92

0,06

24,00

1,02

0,12

27,00

1,04

0,09

30,00

1,00

0,02

33,00

0,97

— 0,01

36,00

0,98

0,00

39,00

0,99

0,02

42,00

1,00

0,02

45,00

1,00

0,00

48,00

1,00

0,00

51,00

1,00

0,00

54,00

1,00

0,00

57,00

1,00

0,00

Таблица 13. Значение переходного процесса в АСР с ПД-регулятором

T

НВ

— 0.97E-02

0.20E-03

0.4

— 0.247 221

— 0.68E-04

0.8

4.27 797

0.30E-02

1.2

11.3696

0.20E-01

1.6

14.6516

0.43E-01

12.7932

0.61E-01

2.4

8.86 159

0.67E-01

2.8

6.40 825

0.65E-01

3.2

6.61 090

0.60E-01

3.6

8.35 361

0.60E-01

9.85 721

0.63E-01

4.4

10.0876

0.67E-01

4.8

9.27 302

0.70E-01

5.2

8.38 501

0.70E-01

5.6

8.10 610

0.68E-01

8.38 829

0.67E-01

6.4

8.83 347

0.67E-01

6.8

9.9 730

0.68E-01

7.2

9.2 298

0.69E-01

7.6

8.75 263

0.69E-01

8.60 058

0.69E-01

Рисунок 17 — Переходный процесс в АСР с ПИ-регулятором по заданию

Рисунок 18 — Переходный процесс в АСР с ПИ-регулятором по возмущению

Рисунок 19 — Переходный процесс в АСР с ПД-регулятором по заданию

Рисунок 20 — Переходный процесс в АСР с ПД-регулятором по возмущению

11. Анализ качества регулирования

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

С помощью переходных кривых (рис. 17−20) определим прямые критерии качества и выберем наилучший закон регулирования.

11.1 Статическая ошибка еуст

Для обоих регуляторов, как по заданию, так и по возмущению еуст > 0.

11.2 Время регулирования

.

11.2.1. Для кривой по задающему воздействию

ПИ-регулятор:; Тр = 36 мин.

ПД-регулятор:; Тр = 60 мин.

11.2.2 Для кривой по возмущающему воздействию

ПИ-регулятор:; Тр = 32 мин.

ПД-регулятор:; Тр = 44 мин.

11.3 Перерегулирование

.

Перерегулирование находится только для процесса по задающему воздействию, так как по возмущению hуст = 0.

ПИ: .

ПД: .

11.4 Степень затухания

: .

11.4.1 Для кривой по задающему воздействию

ПИ: .

ПД: .

11.4.2 Для кривой по возмущающему воздействию

ПИ: .

ПД: .

Анализируя прямые критерии качества можно сделать вывод, что наилучшее регулирование осуществляется в системе с ПИ-регулятором. Поэтому для рассчитанной АСР выбираем ПИ-закон регулирования.

12. Выбор промышленного регулятора

Найдём истинные настройки регулятора. Для этого необходимо учесть коэффициенты усиления датчика и клапана:

;

.

Для нашей АСР выбираем ПИ-регулятор типа ПР 3.31. Его передаточная функция имеет вид:

Найдём значение параметров настройки:

1.Предел пропорциональности:

2.Время изодрома:

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

Список используемой литературы

1. Г. К. Аязян «Расчет автоматических систем с типовыми алгоритмами регулирования» Уфа: УНИ, 1989 г.-135с.

2. Г. К. Аязян «Исследование линейной системы по корневым критериям качества». Методическое руководство. Уфа: УНИ, 1984 г. — 23с.

3. Г. К. Аязян. «Расчет настроечных параметров типовых регуляторов одноконтурных автоматических систем регулирования». Методическое руководство по курсовому и дипломному проектированию. Уфа: УНИ, 1985 г. — 25с.

Показать весь текст
Заполнить форму текущей работой