Лабораторная работа. "Имитационное моделирование системы фазовой автоподстройки частоты в пакете моделирования динамических систем Simulink"

Цель работы: Изучить методы имитационного моделирования системы автоматического регулирования и исследования основных характеристик систем фазовой автоподстройки частоты (ФАП).

Домашнее задание

1.  Составить математическую модель ФАП 1-го порядка и 2-го порядка.

2.  Построить передаточные характеристики ФАП 1-го порядка и 2-го порядка с идеальным и неидеальным интегратором.

3.  Используя передаточные характеристики определить шумовую полосу ФАП 1-го порядка и 2-го порядка и сравнить с расчетными выражениями.

4.  Определить полосу захвата и удержания ФАП 1-го порядка.

Лабораторное задание

Моделирование систем передачи дискретных сообщений (ПДС) и их функциональных узлов на ЭВМ является экономичным способом их исследования и проектирования. Оно позволяет значительно ускорить процесс разработки и найти наиболее оптимальные характеристики их функционирования.

Моделирование на ЭВМ можно условно разбить на следующие этапы:

·  формулировка задачи моделирования, разработка модели системы или устройства и программы их исследования;

·  составление модели;

·  моделирование системы ПДС или устройства при заданных режимах работы и внешних воздействиях;

·  интерпретация результатов моделирования, т.е. построение различных выводов о функционировании объекта по данным, полученным в результате моделирования.

Рис.1. Структурная схема имитационной модели системы ФАП первого, второго и третьего порядка с контрольно-измерительными блоками: Power meter - измеритель мощности; Spectr_a, b, c, d, e - осциллографы и анализаторы спектра

Рис.2 Имитационная модель генератора входного сигнала ФАП с постоянной и линейно-изменяющейся частотной расстройкой.

Рис.3. Имитационная модель фазового детектора.

Рис.4 Имитационная модель ФНЧ фазового детектора.

Рис.5 Имитационная модель фильтра первого и второго порядка кольца ФАП второго и третьего порядка соответственно.

1.  Построить модель фазовой автоподстройки частоты (ФАПЧ) в пакете моделирования динамических систем Simulink. Для этого открыть новое окно модели Simulink. В библиотечном модуле в соответствующих разделах найти типовые функциональные блоки структурной схемы модели ФАП и скопировать их в окно модели. Соединить входы и выходы блоков в соответствии с рис.1.

2.  Настроить общие параметры модели. Для этого необходимо выполнить команду Parameters в позиции Simulation главного меню пакета Simulink. Откроется окно параметров модели. В закладке Solver установить время моделирования (Stop time) равным 5 секунд. Выбрать метод изменения независимых переменных с фиксированным шагом (Type: Fixed-step) и метод решения дифференциальных уравнений при моделировании дискретный (discrete (no continuous states)). Установить время дискретизации (Fixed Step Size) для модели, исходя из следующих предпосылок: для сигналов генераторов с частотой 10 Гц на один период колебания должно приходится 40 отсчетов времени.

3.  Настроить параметры каждого функционального блока структурной схемы. Указать, если требуется, в настройках блока время дискретизации (Sample time). Установить частоты опорного и управляемого генераторов равными 10 Гц (если требуется указать значение в рад/с, ввести выражение 2*pi*f, где f - частота в герцах). Установить разность фаз между фазами опорного и управляемого генераторов (j = jг - jо) равной 90^{о (}pi/2), при которой полезная составляющая на выходе фазового детектора равна нулю. Проконтролировать установленные параметры генераторов, сравнивая осциллограммы их выходных сигналов

4.  Установить постоянную времени RC-фильтра, включенного в блок фазового детектора, обеспечивающую подавление удвоенной частоты входного сигнала на его выходе в 10 раз, снять импульсную и передаточную характеристики фильтра, используя спектроанализатор.

5.  Снять передаточную характеристику идеального пропорционально - интегрирующего фильтра (ПИФ) схемы ФАП, определяемой выражением F (p) = (1 + a/p), где а - коэффициент передачи интегрирующей ветви ПИФ, р - оператор Лапласа. Принять значение а в диапазоне [10,30] с шагом 5.

6.  Оборвать петлю обратной связи ФАП. Подключив измеритель средних значений сигнала к выходу фильтра детектора, путем изменения значения фазы входного сигнала в интервале j Í [p, p], снять дискриминационную характеристику U_д = h (j) фазового детектора.

7.  К входу ГУН подключить источник постоянного напряжения и путем изменения его величины в диапазоне [-1,1] с шагом 0.5 измерить крутизну ГУН и построить характеристику управления Df_у = f (U_у).

8.  Подключив генератор пилообразного напряжения к входу генератора входного сигнала, управляемого напряжением, и установив начальную расстройку Df собственных частот генератора входного сигнала и ГУН, обеспечивающую отсутствие захвата, определить полосу захвата ФАП 1-го и 2-го порядка.

9.  Подключив генератор пилообразного напряжения к входу генератора входного сигнала, управляемого напряжением, и установив нулевую начальную расстройку Df собственных частот генератора входного сигнала и ГУН, определить полосу удержания ФАП 1-го и 2-го порядка.

10.  Анализируя сигнал управления на входе ГУН, определить время ввода в синхронизм ФАП 1-го и 2-го порядка при начальных расстройках в диапазоне значений [-1,1] с шагом 1.

11.  Построить зависимость полосы захвата ФАП 1-го и 2-го порядка от коэффициента усиления петли обратной связи ФАП. Коэффициент усиления менять от 10 до 4 с шагом 2.

12.  Подключить к входу ФАП гармонический сигнал с начальной расстройкой в пределах полосы захвата и сигнал с выхода генератора нормального шума. С помощью измерителя средне - квадратичных значений определить дисперсию сигнала управления и величину фазового джиттера сигнала ГУН, а также определить отношение сигнал-шум внутри кольца ФАП и на ее входе.

Рис.6. Временная диаграмма сигнала управления ГУН ФАП первого порядка при линейном изменении частоты входного сигнала, иллюстрирующая режимы биений, захвата, удержания и срыва синхронизма: Sу = 1 Гц/В, время анализа 40 с, скорость изменения частоты входного сигнала 0.2 Гц/с, полоса удержания равна полосе захвата Δf_уд = Δf_з = 1 Гц, режим слежения наблюдается на временном интервале tÎ [23, 28] с.

Основные определения

Рис. П.1.1 Структурная схема системы фазовой автоподстройки частоты

Рис. П.1.2 Характеристика управления ГУН

Структурная схема системы ФАПЧ представлена на рис.П.1.1 Она состоит из фазового детектора (ФД), фильтра низкой частоты петли регулирования, управляемой системы (УС), включающей генератор, управляемый по частоте напряжением (ГУН), и управляющий элемент (УЭ).

Управляющий элемент предназначен для изменения частоты управляемого генератора.

Зависимость частоты генератора f_г от уровня управляющего напряжения U_у называют характеристикой управления f_г (U_у) (рис.П.1.2). Обычно используют линейный участок этой характеристики. Крутизна управления, или по-другому крутизна управителя, определяется по формуле S_у = Df_г / DU_у и имеет размерность Гц/В.

Фазовый детектор (ФД) создает напряжение, определяемое разностью фаз j колебаний управляемого генератора j_г и колебаний опорного генератора j_о.

Зависимость напряжения U_д на выходе ФД от разности фаз сравниваемых сигналов j = j_г - j_о называют дискриминационной характеристикой фазового детектора U_д (j), а S_д = dU_д/dj - крутизной дискриминационной характеристики, или по-другому крутизной дискриминатора, с размерностью В/Гц. Удобно пользоваться нормированной дискриминационной характеристикой

F (j) = U_д (j) /E,

где E - максимальное напряжение на выходе ФД, которое зависит от амплитуд сравниваемых сигналов и схемы ФД. F (j) - периодическая функция, изменяющая свое значение в пределах ±1.

Характеристика ФД может иметь синусоидальную, рис.П.1.3а, треугольную, рис.П.1.3б, а также знаковую, трапецеидальную или пилообразную формы в зависимости от схемы ФД.

Рис. П.1.3 Нормированные характеристики ФД

 

Фильтр нижней частоты (ФНЧ) применяется для выделения полезной низкочастотной составляющей из сигнала на выходе фазового детектора. В петле ФАП для повышения порядка могут быть применены идеальный и неидеальный интегрирующие фильтры. Неидеально-интегрирующий фильтр характеризуется коэффициентом передачи K (p) = 1/ (1+pT), где T = R*C, рис.1.4а. Пропорционально-интегрирующий фильтр с неидеальным интегратором, рис.1.4б, имеет коэффициент передачи K (p) = (1 + p m T) / (1 + p T), где T = C* (R+R1), m = R1/ (R+R1). Коэффициент передачи пропорционально-интегрирующего фильтра с идеальным интегратором описывается выражением K (p) = ( (а + p) / p), где а - коэффициент передачи интегрирующей ветви.

Рис. П.1.4 Схемы неидеальных ФНЧ:

а) неидеальный интегрирующий фильтр;

б) неидеальный пропорционально-интегрирующий фильтр

Функционирование системы ФАП состоит в следующем. Гармонический сигнал с частотой w_г и фазой j_г с выхода управляемого генератора ФАП поступает на первый вход фазового детектора, на второй вход которого поступает входной гармонический сигнал от опорного генератора с частотой w_о и фазой j_о. Разность частот (f_{г -} f_о) = Δf_н в начальный момент времени определяет начальную частотную расстройку Δf_н системы ФАП. Амплитуда напряжения U_д на выходе ФД зависит от разности фаз сигналов на его входах. Полученное напряжение U_д, после низкочастотной фильтрации, поступает на управляющий элемент, который так изменяет частоту управляемого генератора ФАП, что бы его частота w_г и частота w_о стали равными, а разность фаз на выходе ФД постоянной. Для обеспечения режима захвата и последующего удержания крутизна дискриминатора S_д и крутизна управителя S_у должны быть противоположными по знаку.

Соответственно ФАП, как система автоматического регулирования, является адаптивным (следящим) фильтром, обеспечивающим фильтрацию входного сигнала путем изменения частоты ГУН сигналом управления, формируемого, в общем случае, на выходе фильтра петли ФАП, который подключен к выходу фазового детектора. Для ФАП 1 - го порядка фильтр петли ФАП отсутствует и U_у = U_д. В зависимости от параметров кольца система ФАП может работать как при постоянной так и при изменяющейся частотной расстройке между входным сигналом ФАП и сигналом ГУН.

В результате моделирования должны быть получены качественные и количественные характеристики системы ФАП: полоса захвата, полоса удержания, время вхождения в синхронизм, а также переходные процессы фазы и частоты в системе.

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

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

Если в системе происходит непрерывное нарастание разности фаз управляемого и опорного генераторов, то говорят, что система находится в режиме биений. Переходное состояние системы, при котором режим биений переходит в режим удержания или квазисинхронизма, называется режимом захвата.

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

На практике полоса удержания Δf_уд определяется путем медленного изменения в режиме синхронизма частоты опорного генератора как в сторону повышения, так и в сторону понижения, до моментов выхода системы ФАП из синхронизма. Разность найденных частот характеризует полосу удержания. Соответственно полоса захвата Δf_з определяется по наступлению синхронизма системы ФАП при медленном изменении начальной расстройки Δf_н от больших значений к малым путем изменения частоты опорного генератора. Для ФАП 1 - го порядка Δf_з = Δf_уд = S_уE, причем для шумовых возмущений со стороны входа системы - линеризированная модель ФАП является эквивалентом RC - фильтру с постоянной времени t =1/ S_дS_уE, а установившаяся фазовая ошибка в режиме слежения j_уст = arcsin (Δf_н/ S_уE). Для систем ФАП большего порядка полосы удержания и захвата не равны друг другу.

Элементы теории систем фазового регулирования. Система фазовой автоподстройки частоты

 

Система фазовой автоподстройки частоты состоит из трех основных элементов: перемножителя, инвариантного во времени линейного фильтра и управляемого генератора, рис. П.2.1.

Представим входной сигнал ФАП в виде, соответственно сигнал управляемого генератора ФАП - , где А и K_{1 -} соответственно эффективные напряжения входного сигнала и сигнала генератора управляемого напряжением (ГУН), θ (t) - полная фаза входного сигнала, а θ' (t) - полная фаза сигнала ГУН, причем θ (t) = w₀ t + θ₁ (t) и θ' (t) = w₀ t + θ₂ (t), а θ₁ (t) и θ₂ (t) - составляющие полной фазы относительно линейной фазы.

Еслu отключить управляющий сигнал е (t) от управляемого генератора, то генератор даст синусоидальный сигнал постоянной частоты w₀. После подключения управляющего сигнала частота генератора становится равной [w₀ + K₂e (t)], где K_{2 -} коэффициент пропорциональности размерности радиан на секунду на вольт.

Рис. П.2.1 Система фазовой автоподстройки частоты

Таким образом, частота сигнала ГУН, определяемая производной от полной фазы сигнала на выходе управляемого генератора будет равна

Сигнал x (t) на выходе перемножителя фазового детектора, равный произведению сигналов входного сигнала и сигнала ГУН, определится выражением

Слагаемое суммарной частоты ослабляется совместным действием фильтра и управляемого генератора, так что его можно отбросить.

После прохождения сигналом фазового детектора линейного инвариантного во времени фильтра сигнал е (t) на выходе равен

Причем предполагается, что входной сигнал включен в момент t = 0. Слагаемое e₀ (t) представляет сигнал на выходе, зависящий только от начальных условий в схеме фильтра в момент t = 0. В случае, когда начальные условия можно выбирать произвольно, e₀ (t) º0 при любых значениях t.

Функция h (t) называется импульсной переходной функцией фильтра и представляет обратное преобразование Лапласа от его передаточной функции F (р). В большинстве случаев линейный фильтр представляет электрическую схему, состоящую из линейных элементов с сосредоточенными постоянными (сопротивления, конденсаторы и индуктивности) и, возможно, из линейных усилителей. Соотношение между сигналами на входе и на выходе в такой системе описывается дифференциальным уравнением

где т £ п. Отношение двух многочленов

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

С учетом сделанных обозначений частота сигнала ГУН определится выражением

Определяя фазовую ошибку в виде  а коэффициент усиления петли регулирования ФАП как получим интегро-дифференциальное уравнение, описывающее работу кольца фазовой автоподстройки при известной фазе входного сигнала q (t):

Принимая во внимание, что

уравнение представляется в виде

которому соответствует модель, изображенная на рис.П.2.2.

Рис. П.2.2 Блок-схема системы фазовой автоподстройки частот

В схеме, полученной на основе математической модели, по отношению к схеме на основе физической модели, перемножитель заменяется вычитающим устройством и синусоидальной нелинейностью, а управляемый генератор - интегратором. Последняя из упомянутых замен означает, что фаза выходного сигнала управляемого генератора пропорциональна интегралу управляющего сигнала. Следует заметить также, что усиление петли регулирования возрастает при возрастании амплитуды принятого сигнала. Если фаза принятого сигнала q (t) известна, то фазовую ошибку j (t) можно найти, решив нелинейное интегро-дифференциальное уравнение.

 Линейное приближение и переходный процесс

Если фазовая ошибка j (t) равна нулю, то говорят, что произошел "захват" фазы в системе фазовой автоподстройки частоты.

Если ошибка j (t) все время мала по сравнению с 1 рад, то можно воспользоваться приближением , которое дает ошибку менее 5%. В этом случае говорят, что система регулирования близка к захвату фазы, а синусоидальную нелинейность может быть аппроксимирована линейной зависимостью.

В этом случае работа системы описывается линейным уравнением, которое получается при замене sinφ на φ в интегро-дифференциальном уравнении

которое в предположении, что преобразования Лапласа существуют, в операторной форме имеет вид

где F (p) - передаточная функция линейного фильтра, а j (р) - преобразование Лапласа для j (t), q (p) - преобразование Лапласа для q (t).

Уравнению соответствует блок-схема, изображенная на рис. П.2.3.

j(p)

q1(p)

q2(p)

1/p

F(p)

Рис. П.2.3 Линейная модель системы фазовой автоподстройки частот

Из интегро-дифференциального уравнения в операторной форме нетрудно получить выражения для величины фазовой ошибки φ (р) и фазы θ₂ (р) выходного сигнала ГУН, обеспечивающего слежения за изменениями фазы входного сигнала θ₁ (р)

Отношение Н (р) = θ₂ (р) / θ₁ (р) и равное

называется передаточной функцией замкнутой петли регулирования. Используя ее, получим следующие соотношения между j (p), q₁ (p), q₂ (p):

,

.

Результат обратного преобразования:

называется импульсной переходной функцией замкнутой петли регулирования. Если F (p) является рациональной функцией, то Н (p) - также рациональная функция. Условие устойчивости приводит к требованию, чтобы нули функции [1 + A K F (p) /p] находились в левой полуплоскости. Выполнив обратное преобразование над θ₂ (р) и φ (р), получим

Полученные выражения позволяют исследовать работу ФАП при различных законах изменения фазы входного сигнала θ₁ (р).

Пусть, например, принимаемый сигнал имеет постоянную частоту w [рад/сек] и начальную фазу q₀ и пусть в системе регулирования не будет фильтра. Тогда

.

В этом случае q₂ (t) = (w-w₀) t+q_{0, т}ак что

Обратное преобразование дает

Если предел limj (t) существует, его называют установившейся фазовой ошибкой. Из формулы для j (t) следует, что в рассматриваемом случае установившаяся фазовая ошибка равна (w-w₀) / (AK) [рад], что означает, что управляемый генератор синхронизирован с принимаемым сигналом по частоте, но захват по фазе не может быть достигнут. Для того чтобы линейная модель была применима, необходимо, чтобы величины (w-w₀) / (AK) и q₀ были малы.

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

а преобразование Лапласа для фазовой ошибки соответственно равно

Для определения величины установившейся фазовой ошибки можно воспользоваться предельной теоремой для преобразования Лапласа:

Таким образом, введя в петлю регулирования второй интегратор, можно свести к нулю установившуюся фазовую ошибку в случае принимаемого сигнала, имеющего вид синусоиды постоянной частоты. Захваченная по фазе петля регулирования без фильтра называется петлей регулирования первого порядка, а петля с фильтром, содержащим идеальный интегратор, называется петлей регулирования второго порядка. Вообще, порядок системы регулирования равен числу конечных полюсов передаточной функции разомкнутой системы, т.е. в данном случае числу полюсов функции АКF (p) /p.

Рис. П.2.4 Схема фильтра для ФАП-2

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

Рис. П.2.5 Схема фильтра неидеальной петли ФАП-2

Она совпадает с передаточной функцией системы, состоящей из фильтра низких частот и усилителя и изображенной на рис. П.2.5, где а =1/R₂C и e = l/ (R₁ + R₂) С. Эту схему легче реализовать, чем аналоговый интегратор, изображенный на рис. П.2.4, ее передаточная функция будет близка к передаточной функции аналогового интегратора, если сделать R₁ гораздо больше R₂ и скомпенсировать ослабление, применив усилитель с большим усилением.

Передаточная функция замкнутой системы равна

При такой передаточной функции преобразование Лапласа фазовой ошибки равно

и ее установившееся значение будет

что равно установившейся фазовой ошибке петли первого порядка, уменьшенной в e/а раз. Это и является мерой степени приближения рассматриваемой петли к идеальной петле второго порядка.

Наконец, рассмотрим принимаемый сигнал, частота которого линейно изменяется во времени

где R есть скорость изменения частоты в радианах в секунду за секунду. Это соответствует, например, случаю приема сигнала, передаваемого с помощью генератора постоянной частоты с борта самолета, перемещающегося с постоянным радиальным ускорением Rc/w [м/сек²] по отношению к приемнику, где с - скорость распространения в метрах в секунду. Тогда

Если для слежения за таким сигналом применить петлю первого порядка [F (s) = I], то преобразование Лапласа ошибки будет иметь вид

и фазовая ошибка j (t) неограниченно возрастает при t®¥, как это следует из применения предельной теоремы. Такой же результат получается для неидеальной петли второго порядка. Таким образом, необходимо применить по меньшей мере идеальную петлю второго порядка. В этом случае

так что установившаяся фазовая ошибка будет равна

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

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

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

и преобразование Лапласа фазовой ошибки будет равно

откуда следует, что установившаяся фазовая ошибка равна нулю.

Рассмотренные выше случаи сведены в табл.1. Из нее видно, в частности, что для отслеживания постоянной частоты (или угловой скорости) с конечной ошибкой достаточно применить петлю первого порядка, а для отслеживания линейно изменяющейся частоты (или углового ускорения) с конечной ошибкой необходима петля второго порядка. Увеличение на единицу порядка системы приводит к устранению установившейся ошибки, а понижение порядка на единицу приводит к неограниченному возрастанию ошибки.

Эти замечания справедливы для всех линейных систем регулирования.

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

Рис. П.2.6 Схема фильтра петли регулирования ФАП-3

Таблица 1

Фаза принимаемого сигнала	Порядок системы	Передаточная функция фильтра	Передаточная функция замкнутой системы	Установившаяся ошибка	Шумовая полоса ФАП
wt+q0	первый	1	AK/ (p+AK)	(w-w0) /AK	AK/4
wt+q0	второй	1+a/p	AK (p+a) / (p2 +Akp+Aka)	0	(AK+а) /4
wt+q0	второй (неидеальной)	(p+a) / (p+e)	AK (p+a) / (p2 + (Ak+e) +Aka)	(e/a) ( (w-w0) /AK)	(АК/4) * [ (AK+а) / (AK+ε)]
1/2 (Rt2+wt+q0)	второй	1+a/p	AK (p+a) / (p2 +Akp+Aka)	R/aAK	(AK+а) /4
1/2 (Rt2+wt+q0)	Третий	1+a/s+b/ s2	AK (s2+as+b) / (s3 +Aks2+aAks+bAK)	0	(АК/4) * [ (аAK+а2-b) / (aAK-b)]

Литература

1.  Радиотехнические цепи и сигналы. Под ред. К.А. Самойло. - М.: Радио и связь, 1982.

2.  Э.Д. Витерби. Принципы когерентной связи. - М.: Советское радио, 1970.

3.  Г.А. Емельянов, В.О. Шварцман. Передача дискретной информации. - М.: Связь, 1982.

4.  В.С. Чернега. Проектирование технических средств обмена и передачи информации. - М.: Высшая школа, 1990.

Д.А. Абдулаев, М.Н. Арипов. Передача дискретных сообщений в задачах и упражнениях. - М.: Радио и связь, 1985.

Лабораторная работа. "Имитационное моделирование системы фазовой автоподстройки частоты в пакете моделирования динамических систем Simulink" Цель работы: Изучить методы имитационного моделирования системы автоматического регулирова

 

• 14:21 05 Июля 2016
  Российско-китайский медуниверситет откроет магистерскую программу по TCM
• 01:21 05 Июля 2016
  На форуме МГУ Россия и Китай обсудят научно-образовательное сотрудничество
• 11:00 04 Июля 2016
  Реморенко: "ЕГЭ-туризм" в России сходит на нет
• 10:52 04 Июля 2016
  Летняя российско-австрийская школа откроется в Москве
• 15:31 01 Июля 2016
  В СПбГУ могут воссоздать географический факультет

Заказать уникальную работу