Модель Леонтьева многоотраслевой экономики

 

1. Анализ объекта управления. Постановка задачи

макроэкономический модель контур оценивание

Рассматривается макроэкономическая модель многоотраслевой экономики В.В. Леонтьева.

, (количество отраслей, )

- мощность производительных фондов

- инвестиции в отрасль

- коэффициенты амортизационных отчислений

Балансные соотношения (учитывают взаимосвязи между отраслями):

- суммарные инвестиции

- коэффициент распределения инвестиций

- коэффициенты, характеризующие взаимосвязи между отраслями

Эффективные производственные мощности:

- коэффициенты использования

- суммарная мощность производства

Суммарные продукт:

Области изменения коэффициентов:

(разнотемповые процессы, отличие на порядок)

Структурная схема:

Схема

Выберем коэффициенты :

2. Математические модели объекта управления

.1 Уравнения в переменных состояния

Система линейных дифференциальных уравнений, которые описывают систему:

Уравнения в переменных состояния:

Тогда запишем компоненты этих уравнений:

2.2 Операторное описание через передаточную функцию

Найдём передаточную функцию системы, устанавливающую связь между входом и выходом, используя предыдущее описание:

Здесь - союзная матрица (матрица, состоящая из алгебраических дополнений).

Вычисляя передаточную функцию системы с помощью MatLab (см. файл №1 в приложении), получаем тот же самый результат:

H =(-k*a2*k2+k*k1*p+k*k1*mu2+k2*p-k2*p*k+k2*mu1-k2*mu1*k-a1*k1+a1*k1*k)/p/(p^2+p*mu2+mu1*p+mu1*mu2-a1*a2)

Найдём передаточные функции каждой из двух отраслей:

Применяя преобразования Лапласа к этим уравнениям, получим:

Находя отсюда и , получаем передаточные функции и для каждой из отраслей и :

2.3 Уравнения вход-выход

Используя передаточную функцию, полученную в предыдущем параграфе, опишем систему в уравнениях вход-выход:

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

2.4 Описание в виде интеграла свёртки

Найдём весовую функцию . Так как рассматриваемая система является стационарной, то эта весовая функция будет иметь один аргумент и определяться через обратное преобразование Лапласа (или матричную экспоненту):

Для упрощения вида этой функции введём обозначения - корень числителя и два ненулевых корня знаменателя соответственно, взятые с обратными знаками:

Разложим передаточную функцию на простые дроби:

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

Также не будем здесь приводить результат, полученный с помощью MatLab с использованием символьной алгебры и матричной экспоненты (см. файл №2 в приложении), так как он тоже очень громоздок и неинформативен.

2.5 Частотные характеристики

Для определения частотных характеристик системы достаточно подставить в передаточную функцию аргумент вместо :

Найдём частотные характеристики с помощью MatLab (см. файл №4 в приложении) и убедимся, что результаты совпадают:

H =-(-k*a2*k2+i*k*k1*w+k*k1*mu2+i*k2*w-i*k2*w*k+k2*mu1-k2*mu1*k-a1*k1+a1*k1*k)/w/(i*w^2+w*mu2+w*mu1-i*mu1*mu2+i*a2*a1)

Найдем модуль и аргумент этой функции:

3. Свойства объекта управления

.1 Устойчивость

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

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

Рассмотрим систему без интегратора:

Рис.

Знаменатель передаточной функции

Для полинома второй степени необходимый и достаточный критерий устойчивости Гурвица эквивалентен необходимому условию устойчивости Стодолы:

Так как существует ограничение на коэффициенты (они должны быть положительными), то второе условие системы выполняется. Остается последнее условие:

Построим корневой годограф с помощью MatLab (см. файл №4 в приложении), стрелками покажем направление движения корней характеристического полинома с ростом величины произведения :

Рис.

3.2 Анализ минимальнофазовости объекта

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

Так как функции MatLab, служащие для определения свойств объекта управления, не работают с символьными переменными, то будем использовать числовые значения параметров, которые будут получены в §3.5. Исследуем минимальнофазовость (см. файл №5 в приложении):

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

3.3 Исследование управляемости и наблюдаемости

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

Потребовав, чтобы корни числителей не равнялись корням знаменателя, получим серию ограничений на коэффициенты системы:

Так как функции MatLab, служащие для определения свойств объекта управления, не работают с символьными переменными, то будем использовать числовые значения параметров, которые будут получены в §3.5. Исследуем управляемость и наблюдаемость (см. файл №6 в приложении):

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

3.4 Анализ установившихся режимов

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

3.5 Окончательный выбор параметров и его обоснование

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

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

Теперь подставим их в систему ограничений (1)-(12):

Отобразим все эти условия на плоскости с помощью MatLab (см. файл №7 в приложении), причем в условиях (2) - (11) заменим знаки и на равенство:

Рис.

На данном графике изображены: условие (2) - розовым цветом, условие (3) - зелёным цветом, условие (4) - синим цветом, условия (5) - (11) - чёрным цветом и условие (12) - красной пунктирной линией.

Для большей информативности рассмотрим этот график при изменениях параметров , :

Рис.

Исходя из этого графика выберем параметры и .

Условие (1) - положительность коэффициентов и - учтено на графике, так как мы и рассматриваем его при положительных параметрах.

Условие (2) эквивалентно тому, что область изменения коэффициентов и лежит ниже розовой кривой.

Условие (3) эквивалентно тому, что область изменения коэффициентов и лежит левее зелёной прямой.

Условие (4) эквивалентно тому, что область изменения коэффициентов и лежит ниже синей прямой.

Условия (5) - (11) эквивалентны тому, что область изменения коэффициентов и не должна содержать чёрные кривые с неким запасом, т.е. окрестностью. Другими словами, коэффициенты и не должны лежать на чёрных кривых или вблизи них.

Условие (12) эквивалентно тому, что коэффициенты и должны лежать на красной пунктирной прямой.

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

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

Итак, запишем все выбранные параметры системы:

4. Процессы в объекте управления

.1 Процессы при импульсном воздействии

Рассмотрим для начала два инерционных звена нашей системы:

Быстрота протекания процессов и время установления:

Графики реакции на импульсное воздействие инерционных звеньев (см. файл №8 в приложении):

Рис.

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

Теперь рассмотрим реакцию на импульсное воздействие как всей системы, так и каждой из её отраслей (см. файл №8 в приложении):

Рис.

Рис.

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

На втором графике установившийся уровень ненулевой, что объясняется наличием интегратора в системе.

Вычислим значения реакции в начальный момент времени и установившиеся уровни и убедимся в том, что они соответствуют графикам:

4.2 Процессы при ступенчатом воздействии

Рассмотрим два инерционных звена нашей системы:

Построим графики реакции этих звеньев на ступенчатое воздействие (см. файл №8 в приложении) и рассчитаем их установившийся уровень:

Рис.

Установившийся уровень:

Тогда

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

Теперь рассмотрим реакцию отраслей системы на ступенчатое воздействие (см. файл №8 в приложении):

Рис.

Установившийся уровень для двух отраслей один и тот же (мы сами потребовали это в §3.4):

Теперь рассмотрим реакцию всей системы на ступенчатое воздействие (см. файл №8 в приложении):

Рис.

Этот график не ограничен:

Теперь рассчитаем некоторые показатели для переходных процессов:

Показатели, характеризующие колебательность системы - нулевые, т.к. корни знаменателя передаточной функции - вещественные:

Показатель колебательности:

Склонность системы колебаниям:

Показатель степени устойчивости:

Время установления: (время вхождения в пятипроцентную область, )

Мера быстродействия:

Стоит добавить, что система обладает астатизмом первого порядка, так как она имеет один интегратор.

Приведём здесь корневой годограф системы с окончательно определёнными параметрами (см. файл №8 в приложении):

Рис.

4.3 Процессы при гармоническом воздействии

Сначала рассмотрим графики логарифмических амплитудно- и фазово-частотных характеристик инерционных звеньев (см. файл №8 в приложении):

Рис.

Установившаяся реакция системы на гармоническое воздействие будет тоже гармонической, но с другой амплитедой и сдвигом по фазе:

Построим графики логарифмических амплитудно- и фазово-частотных характеристик и годограф (график полной частотной характеристики) двух подсистем и системы в целом (см. файл №8 в приложении):

Рис.

Рис.

Рис.

Рис.

Так как система и её составляющие физически реализуемы, то АЧХ должна убывать, ФЧХ должна быть отрицательной (реакция всегда запаздывает по отношению к воздействию), а годограф должен заканчиваться в точке . Все эти требования как раз и отражены на приведённых выше графиках.

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

.1 Выбор метода и построение структурной схемы с обратной связью

Будем строить закон управления системой по состоянию с оптимизацией контура управления:

Схема

5.2 Оптимизация контура управления

Мы хотим построить управление, зависящее от состояния так, чтобы замкнутая система была оптимальной в смысле некоторого показателя качества:

,

Где - оценка вектора состояния, - управление, зависящее от желаемого выхода, а - искомая матрица (в нашем случае - строка) коэффициентов обратной связи.

Цель управления: отработка импульсных инвестиций с минимальными затратами на управление.

Итак, интегрально-квадратичный показатель качества:

Так как пара матриц - управляемая, то существует симметричная положительно определённая матрица и матрица коэффициентов обратной связи такие, что при любых начальных условиях оптимальная стабилизация системы обеспечивается линейным управлением по состоянию, которое даёт минимальное значение показателя качества , причём матрица находится решением алгебраического уравнения Лурье-Риккати:

Найдём это управление с помощью встроенных функций MatLab (см. файл №9 в приложении), при этом меняя элементы матриц и будем анализировать собственные числа матрицы замкнутой системы . Таким образом, варьируя параметры и , выберем желаемые собственные числа с небольшой мнимой частью (чтобы в системе была малая колебательность) и по возможности лежащие левей собственных чисел объекта управления (чтобы управление успевало сформироваться, пока в объекте протекают процессы):

= 0.0500= 1.0e-005 *

0.0100 0 0

0.0100 0

0 1.0000

K = 0.1813 -0.0687 0.0141

P =

.1444 -0.1939 0.0014

.1939 0.2822 -0.0004

.0014 -0.0004 0.0002

e =

.0019

.0947 + 0.0565i

.0947 - 0.0565i

Таким образом, полученная строка обратных связей, матрица замкнутой системы и её собственные числа:

Рис.

Показатель степени устойчивости (косвенно характеризует быстродействие):

- корень, лежащий ближе других к мнимой оси - даёт самый медленный процесс.

5.3. Настройка контура оценивания

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

( - матрица невязки)

(ошибка оценивания)

Вычтем из уравнения объекта уравнение оценки:

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

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

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

В качестве желаемого расположения собственных чисел матрицы возьмём распределение Баттерворта:

Корни такого полинома размещаются в вершинах правильного -угольника, а число определяет радиус их распределения.

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

Характеристическое уравнение матрицы :

Приравнивая коэффициенты получившегося полинома к коэффициентам многочлена Баттерворта, получаем систему уравнений для компонентов столбца :

Отсюда решение:

Проверим получившийся ответ: пусть , тогда:

Характеристическое уравнение матрицы :

,

что и соответствует полиному Баттерворта третьей степени (с небольшой погрешностью).

Итак, строим нашу систему и смотрим, как фильтруется возмущение:

Заметим, что получившаяся система имеет порядок .

Итак, начинаем подбирать , на вход системы подаём воздействие с возмущением в виде случайного процесса (белый шум с нулевым математическим ожиданием и заданным отношением сигнал/шум) (см. файл №10 в приложении):

В качестве основного входного воздействия возьмём гармоническое воздействие:

Подаём на вход гармоническое воздействие без помех:

Рис.

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

Рис.

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

Рис.

Подаём на вход гармоническое воздействие без помех:

Рис.

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

Рис.

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

Рис.

Подаём на вход гармоническое воздействие без помех:

Рис.

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

Рис.

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

Рис.

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

Итак, окончательный выбор параметров и завершение построения системы:

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

Собственные числа матрицы :

Незначительные отклонения от желаемых собственных чисел вызваны округлениями в процессе расчётов.

Корневой годограф:

Рис.

5.4 Окончательные результаты и вывод

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

Для каждой из отраслей:

)без обратной связи:

Рис.

)с обратной связью:

Рис.

Для всей системы:

1)без обратной связи:

Рис.

2)с обратной связью:

Рис.

Как мы видим, для каждой из отраслей введение контура управления и оценивания практически не повлияло на время установления, зато добавило в систему колебательность (из-за мнимых частей собственных чисел). Это обусловлено тем, что контур управления строился с учётом ограничения на управление (вкладываемые средства). Если бы мы задавали этот контур по желаемым собственным числам, то можно было бы значительно улучшить процессы в отраслях (т.е. убрать колебательность и уменьшить время установления).

Для всей же системы процессы стали намного лучше: время установления уменьшилось на порядок (), а установившийся уровень стал нулевым.

Вывод

Задав систему (макроэкономическую модель многоотраслевой экономики В.В. Леонтьева), мы исследовали её различные свойства и исходя из них выбрали её некоторые параметры. В дальнейшем мы анализировали процессы в системе при типовых воздействиях. Поставив цель - стабильная работа системы при ограниченном ресурсе - мы построили закон управления системой, оптимизировали контур управления и оценивания. Также система была промоделирована на ЭВМ (с помощью MatLab). В результате мы улучшили процессы в системе при однократном вложении в неё средств (импульсном воздействии). Таким образом можно сравнить, как будет работать система при управлении человеком напрямую и при присутствии в ней внутренних обратных связей.

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

.Курс теории автоматического управления: учеб. пособ. - Первозванский А.А., М.: Наука. Гл. ред. физ.-мат. лит., 1986

.Теория автоматического управления. Линейные системы. - Мирошник И.В., СПб.: Питер, 2005

.Теория автоматического управления в примерах и задачах с применением пакета Matlab. Учебное пособие. - Л.В.Бабко, В.П.Васильев,

.В.С.Королев, Н.Д.Тихонов , Санкт-Петербургский государственный технический университет, СПб., 2001.

Приложение

Файл №1 (вычисление передаточной функции системы)

syms mu1 mu2 a1 a2 k1 k2 k p;=[-mu1 -a1 0; -a2 -mu2 0; k1 k2 0];=[k; 1-k; 0];=[0; 0; 1];=eye(3);=(C')*(p*E-A)^(-1)*B;=simplify(HH)

Файл №2 (вычисление весовой функции)mu1 mu2 a1 a2 k1 k2 k p t;=[-mu1 -a1 0; -a2 -mu2 0; k1 k2 0];=[k; 1-k; 0];=[0; 0; 1];=(C')*expm(A*t)*B

Файл №3 (определение частотных характеристик)

syms mu1 mu2 a1 a2 k1 k2 k w;

A=[-mu1 -a1 0; -a2 -mu2 0; k1 k2 0];=[k; 1-k; 0];=[0; 0; 1];=eye(3);=(C')*(i*w*E-A)^(-1)*B;

H=simplify(HH)

Файл №4 (построение корневого годографа)

a=-1:0.005:1;i=1:401(i)=(-0.11-sqrt(0.0081+4*a(i)))/2;(i)=(-0.11+sqrt(0.0081+4*a(i)))/2;(real(p1(i)),imag(p1(i)), 'ro',real(p2(i)),imag(p2(i)), 'bo');

xlabel('Re(p)');('Im(p)');

hold on

Файл №5 (исследование минимальнофазовости объекта)

A=[-0.01 -0.1478 0; -0.005 -0.1 0; 0.9 0.8 0];=[0.6; 0.4; 0];=[0; 0; 1]';=0;=ss(A,B,C,D);=zpk(sys)

Файл №6 (исследование управляемости и наблюдаемости объекта)

A=[-0.01 -0.1478 0; -0.005 -0.1 0; 0.9 0.8 0];=[0.6; 0.4; 0];=[0; 0; 1]';=0;=ss(A,B,C,D);=ctrb(sys)=obsv(sys)(Co)(Oo)

Файл №7 (построение условий для коэффициентов перекрёстных связей)

s2=inline('a1*a2-0.001','a1','a2');=inline('a1-0.15','a1','a2');=inline('a2-0.0067','a1','a2');=inline('(0.667*a1+sqrt(a1*a2+0.0002025)-0.045)','a1','a2');=inline('(0.667*a2-sqrt(a1*a2+0.0002025)-0.045)','a1','a2');=inline('(1.5*a2+sqrt(a1*a2+0.0002025)+0.045)','a1','a2');=inline('(1.5*a2-sqrt(a1*a2+0.0002025)+0.045)','a1','a2');=inline('(0.558*a2+0.419*a1-0.0665)','a1','a2');=inline('(0.558*a2+0.419*a1-0.01+sqrt(a1*a2+0.0002025))','a1','a2');=inline('(0.558*a2+0.419*a1-0.01-sqrt(a1*a2+0.0002025))','a1','a2');=inline('0.48*a2-0.36*a1+0.0508','a1','a2');on;(s2, 'm-','LineWidth',2);(s3, 'g-','LineWidth',2);(s4, 'b-','LineWidth',2);(s5, 'k-','LineWidth',2);(s6, 'k-','LineWidth',2);(s7, 'k-','LineWidth',2);(s8, 'k-','LineWidth',2);(s9, 'k-','LineWidth',2);(s10, 'k-','LineWidth',2);(s11, 'k-','LineWidth',2);(s12, 'r--','LineWidth',2);

axis equal;

Примечание: в данном файле используется функция implot2.m, служащая для построения графиков функций, заданных в неявном виде; данная функция доступна на официальном сайте MatLab.

Файл №8 (импульсное, ступенчатое, и гармоническое воздействия)

in1=tf([1],[1 0.01])=tf([1],[1 0.1])(in1, 'r')onon(in2, 'b')

legend('1/p+0.01','1/p+0.1')

sys1=tf([0.54 0.000792],[1 0.11 0.000261])=tf([0.32 0.0008],[1 0.11 0.000261])(sys1, 'r')onon(sys2, 'b')('H1(p)','H2(p)')=tf([0.86 0.001592],[1 0.11 0.000261 0])(sys, 'b')on('H(p)')=tf([1],[1 0.01])=tf([1],[1 0.1])(in1, 'r')onon(in2, 'b')

legend('1/p+0.01','1/p+0.1')

sys1=tf([0.54 0.000792],[1 0.11 0.000261])=tf([0.32 0.0008],[1 0.11 0.000261])(sys1, 'r')onon(sys2, 'b')('H1(p)','H2(p)')=tf([0.86 0.001592],[1 0.11 0.000261 0])(sys, 'b')on('H(p)')=tf([0.86 0.001592],[1 0.11 0.000261 0])(sys, 'b')on('H(p)')=tf([1],[1 0.01])=tf([1],[1 0.1])(in1, 'r')onon(in2, 'b')

legend('1/p+0.01','1/p+0.1')

sys1=tf([0.54 0.000792],[1 0.11 0.000261])=tf([0.32 0.0008],[1 0.11 0.000261])(sys1, 'r')onon(sys2, 'b')('H1(p)','H2(p)')=tf([0.86 0.001592],[1 0.11 0.000261 0])(sys, 'b')on('H(p)')=tf([0.54 0.000792],[1 0.11 0.000261])=tf([0.32 0.0008],[1 0.11 0.000261])(sys1, 'r')onon(sys2, 'b')('H1(p)','H2(p)')=tf([0.86 0.001592],[1 0.11 0.000261 0])(sys, 'b')on('H(p)')

Файл №9 (решение задачи оптимальной стабилизации)

A=[-0.01 -0.1478 0; -0.005 -0.1 0; 0.9 0.8 0];=[0.6; 0.4; 0];=0.05;=[10^(-7) 0 0;0 10^(-7) 0;0 0 10^(-5)]

[K, P, e]=lqr(A, B, Q, R)

Файл №10=[-0.01 -0.1478 0; -0.005 -0.1 0; 0.9 0.8 0];=[0.6; 0.4; 0];=[0 0 1]';=[0.1813 -0.0687 0.0141];=0.5;=[4.4122*r^3+1.3767*r^2-0.1517*r+0.0083 -4.9637*r^3+0.9512*r^2-0.10435*r+0.00547 2*r-0.11]';=[A-B*K B*K; zeros(3) A-L*(C')];=[B; 0; 0; 0];=[C;0;0;0];=ss(A1,B1,C1',0)

t=0:0.1:100;=sin(0.3*t);=awgn(x,1);

y2=awgn(x,10);(sys, x, t)(sys,y1,t)(sys,y2,t)=tf(sys)(A1)

1. Анализ объекта управления. Постановка задачи макроэкономический модель контур оценивание Рассматривается макроэкономическая модель многоотраслевой экон

Больше работ по теме: