Модель Леонтьева многоотраслевой экономики
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 в приложении):
Рис.
Первый процесс устанавливается значительно
дольше второго, что и определяется временем установления.
Рис.
Рис.
Первый график показывает, что процессы
разнотемповые, но установившийся уровень у них один и тот же.
На втором графике установившийся уровень
ненулевой, что объясняется наличием интегратора в системе.
Вычислим значения реакции в начальный момент
времени и установившиеся уровни и убедимся в том, что они соответствуют
графикам:
4.2 Процессы при ступенчатом
воздействии
Рассмотрим два инерционных звена нашей системы:
Построим графики реакции этих звеньев на
ступенчатое воздействие (см. файл №8 в приложении) и рассчитаем их
установившийся уровень:
Рис.
Установившийся уровень:
Тогда
Итак, в силу разнотемповости процессов мы снова
видим, что время установления первого процесса значительно больше времени
установления второго процесса, а установившиеся уровни как раз равны этому
времени установления для каждого процесса соответственно.
Теперь рассмотрим реакцию отраслей системы на
ступенчатое воздействие (см. файл №8 в приложении):
Рис.
Установившийся уровень для двух отраслей один и
тот же (мы сами потребовали это в §3.4):
Теперь рассмотрим реакцию всей системы на
ступенчатое воздействие (см. файл №8 в приложении):
Рис.
Этот график не ограничен:
Теперь рассчитаем некоторые показатели для
переходных процессов:
Показатели, характеризующие колебательность
системы - нулевые, т.к. корни знаменателя передаточной функции - вещественные:
Показатель колебательности:
Склонность системы колебаниям:
Показатель степени устойчивости:
Время установления: (время
вхождения в пятипроцентную область, )
Мера быстродействия:
Стоит добавить, что система обладает астатизмом
первого порядка, так как она имеет один интегратор.
Приведём здесь корневой годограф системы с
окончательно определёнными параметрами (см. файл №8 в приложении):
Рис.
4.3 Процессы при гармоническом
воздействии
Сначала рассмотрим графики логарифмических
амплитудно- и фазово-частотных характеристик инерционных звеньев (см. файл №8 в
приложении):
Рис.
Установившаяся реакция системы на гармоническое
воздействие будет тоже гармонической, но с другой амплитедой и сдвигом по фазе:
Построим графики логарифмических амплитудно- и
фазово-частотных характеристик и годограф (график полной частотной
характеристики) двух подсистем и системы в целом (см. файл №8 в приложении):
Рис.
Рис.
Рис.
Рис.
Так как система и её составляющие физически
реализуемы, то АЧХ должна убывать, ФЧХ должна быть отрицательной (реакция
всегда запаздывает по отношению к воздействию), а годограф должен заканчиваться
в точке .
Все эти требования как раз и отражены на приведённых выше графиках.
5. Синтез законов управления для
линейных динамических систем
Будем строить закон управления системой по
состоянию с оптимизацией контура управления:
Схема
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
. В.С.Королев,
Н.Д.Тихонов , Санкт-Петербургский государственный технический университет,
СПб., 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)