«Председатель: д. т. н., профессор И. В. Черных Секция 5. Моделирование в Simulink Оглавление Кудашов А. В., Чернецов В. И. МОДЕЛИРОВАНИЕ АЦП ДЛЯ ИЗМЕРЕНИЙ ПОКАЗАТЕЛЕЙ КАЧЕСТВА ЭЛЕКТРОЭНЕРГИИ Курицына В. В. ...»

Есть альтернативная точка зрения, которая утверждает, что то, что внешне выглядит как динамика трения, на самом деле является результатом высокочастотных нормальных перемещений, вызываемых столкновениями неровностей поверхностей. Эту точку зрения высказал ещё в 1967 г.

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

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

2. Трение как статическая нелинейность (статические модели трения) Статические модели трения (СМТ) не имеют внутренних состояний и характеризуются неизменной зависимостью момента трения от скорости.

Рассмотрим их по возрастанию сложности.

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

2.1. Кулоново трение Характеристика Кулонова трения описывается функцией «сигнум» в виде где y —сила F (момент М) трения, x — скорость , k — коэффициент, определяющий уровень Кулонова трения. Далее не будем акцентировать внимание на том, какой вид движения рассматривается — поступательное или вращательное. Поэтому под «скоростью» будем понимать как линейную, так и угловую скорости.

Составим модель данной функции и проверим ее работу (рис. 5).

Блок «Sign» является стандартным блоком системы Simulink. Однако данный блок обеспечивает амплитуду выходного сигнала, равную только единице, что не достаточно для построения требуемой характеристики, например, «кулонова» трения неединичного значения. Поэтому последовательно с ним включим блок k.

На вход модели следует подать сигнал, близкий к характеристике изменения скорости в современных электромеханических системах. Выберем синусоиду с единичной амплитудой и частотой 1 рад/с. Ее вид представлен на рис. 6. Данный сигнал будем использовать в качестве входного воздействия для испытания всех остальных моделей трения.

Рис. 6. Входное воздействие для испытания моделей моментов трения.

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

Рис. 7. Характеристики сухого трения: а) сигналы «результат» и «момент трения»

при k=1; б) сигналы «результат» и «момент трения» при k=2.

На основании моделирования можно заключить, что с помощью выражения (1) может быть задано только Кулоново трение.

2.2. Кулоново + вязкое трение Характеристика рассматриваемой нелинейности описывается следующим уравнением где y — момент трения М, x — скорость, k1 — коэффициент вязкого трения, k2 — уровень Кулонова трения.

Для моделирования характеристики трения «Кулоново + вязкое трение» в Simulink имеется готовый блок «Coulomb & Viscous Friction» (рис.

8), работа которого построена на основе выражения (2).

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

Рис. 8. Схема моделирования «Кулонова + вязкого» трения.

Данный блок имеет два изменяемых параметра, одним из которых — «Coulomb friction value (Offset)» — задается уровень Кулонова трения, а другим — «Coefficient of viscous friction (Gain)» — коэффициент вязкого трения. На рис. 9 изображены результаты моделирования при единичных значениях указанных параметров.

Рис. 9. Характеристики модели «Кулоново + вязкое трение»

Недостатками обеих рассмотренных моделей является следующее.

Линейного нарастания силы трения от нуля не происходит. Сила трения мгновенно принимает значение Кулонова трения. При моделировании сила трения занимает неопределенное положение между положительным и отрицательным уровнями характеристики трения.

2.3. Модели трения с эффектом Штрибека Наиболее полным математическим описанием характеристики трения, представленным в большинстве как отечественных, так и зарубежных работ до начала 90-х годов, являются различные описания кривой Штрибека в виде статической зависимости. Например, существует Лоренцианова (Lorentzian) и Гауссова (Haussian) модели, модель Тастина (Tustin). В отечественной литературе математические описания кривой Штрибека представлены, в частности, в работах [6], [10], [11].

В [10] предлагается следующая зависимость между моментом трения и числом оборотов двигателя:

где m — коэффициент крутизны наклона спадающей части характеристики, С — уровень трения Кулона, а С2 — участок характеристики, характеризующий уровень статического трения и величину участка отрицательного вязкого трения. Они задаются так, как показано на рис. 10.

Рис. 10. Характеристика и параметры модели (3).

На основе формулы (3) построим модель трения (рис. 11).

Модель имеет в своем составе стандартный блок «Switch». Свернём схему на рис. 11 в один блок (рис. 12) и промоделируем её (рис. 13).

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

Рис. 12. Схема моделирования трения на основе формулы (3).

Рис. 13. Характеристики модели трения на основе формулы (3).

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

Рассмотрим модель трения, предложенную в [11]. Здесь предлагается задавать верхнюю часть характеристики соотношением а нижнюю часть характеристики - соотношением Для их реализации составим следующую структурную схему — рис. 14.

Рис. 14. Схема моделирования трения на основе модели из [11].

Объединив для удобства данную схему в один блок, как показано на рис. 15, и промоделировав ее, получим характеристики, изображенные на рис. 16.

Рис. 16. Характеристики модели трения из [11].

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

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

Рассмотрим более подробно реализацию характеристики момента трения с помощью данной модели.

Рис. 17. Модель из [11] в оригинальных обозначениях.

Параметрами FH± задается уровень статического трения, параметрами F± — уровень трения Кулона (FH±-F± — значение трения Кулона), 0± — коэффициенты крутизны спадающей части характеристики (отрицательное вязкое трение) — чем меньше коэффициент, тем круче падение графика, b± — коэффициенты вязкого трения, — чем больше коэффициент, тем резче подъем графика. Величина петли гистерезиса прямо пропорциональна амплитуде поданного сигнала.

В [11] приведены конкретные значения параметров для каждой кривой характеристики. Представим их:

FH+=1,79 H, F+=0,97 H, 0+=0,003915 м/c, b+=0 Hc/м, FH-=-1,63 H, F-=-0,84 H, 0-=-0,003915 м/c, b-=0 Нс/м.

Введем данные значения параметров в модель на рис. 15 и промоделируем ее. Результаты представлены на рис. 18.

Рис. 18. Характеристики модели трения [11] с параметрами из статьи.

Существует ряд статических моделей трения, задачей которых, в общем, является описание кривой Штрибека, а режим предварительного смещения характеризуется либо функцией sgn, либо одной из сглаживающих функций [12]. Наиболее широко принятой зависимостью для описаСекция 5. Моделирование в Simulink ния статической характеристики момента трения является следующая модель в которой FC — уровень Кулонова трения, FS — уровень статического трения при трогании с места, 2 — коэффициент вязкого трения для высоких скоростей, S — скорость Штрибека (крутизна падения момента трения в начале движения). Таким образом, модель характеризуется четырьмя «статическими» параметрами. Это выражение используется также в современных динамических моделях трения.

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

Основными особенностями динамических моделей являются устранение неопределенности в нулевой области, моделирование предварительного смещения, учет переменного статического трения и петли гистерезиса в кривой Штрибека. Обзор динамических моделей трения и их сравнительное исследование выполнены в [12].

Первую динамическую модель трения предложил, судя по статье [9], Даль (Dahl, 1968). В [13] упоминаются еще несколько ДМТ, предложенных Bo и Pavelescu (1982), Armstrong-Helouvry (1991), Harnoy и Friedland (1994), Canudas de Wit, Olsson, Astrm и Lischinsky (1995).

Bo и Pavelescu предлагают модель, в которой есть две новых величины (dxS /dt и ) в показателе экспоненты Штрибека (см. с. 25 [13]):

В зависимости от выбора значений параметров FC, FS, x S и, модель даёт множество кривых различного вида.

3.1. Модель ЛуГре Наиболее широкое распространение получила в последние годы динамическая модель ЛуГре [14] (oт названий городов в которых работают авторы [9] Лунд (K.J. strm) и Греноблъ (C. Canudas de Wit и P.

Lischinsky), в англоязычной литературе LuGre от Lund-Grenoble):

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

где В данной модели: F — сила трения, — скорость, z — переменная, характеризующая отклонение микронеровностей, расположенных между трущимися поверхностями, 0 - коэффициент упругой жесткости, характеризующий пружиноподобное поведение системы с трением для малых смещений, 1 — коэффициент вязкого трения для скоростей, близких к нулю, 2 — коэффициент вязкого трения для высоких скоростей, FC — уровень Кулонова трения, FS — уровень статического трения, S — скорость Штрибека. Таким образом, данная модель характеризуется шестью параметрами.

На основе выражений (4)–(6) построим модель, структурная схема которой представлена на рис. 19.

Рис. 19. Структурная схема динамической модели трения.

Структурная схема блока 1/g(v) изображена на рис. 20.

Используем для моделирования значения параметров, представленные в [9]. В [9] утверждается, что модель, основывающаяся на выражениях (4) — (6), включает в себя все динамические свойства трения, такие как смещение перед скольжением, фрикционное запаздывание, переменная сила трения при отрыве и «стик-слип» эффект. Исследуем эти свойства с помощью моделирования.

Для демонстрации смещения перед скольжением в [9] проводится следующий опыт. К телу единичной массы была приложена внешняя сила, которая была медленно увеличена до значения 1.425 Н, что составляет 95% от FS. Некоторое время эта сила поддерживалась равной значению 1.425 Н, затем снижена до величины —1.425 Н, некоторое время поддерживалась равной этому значению и затем была снова увеличена до значения 1.425 Н.

Результат представлен в виде графика зависимости силы трения от смещения. Выполним данный опыт с помощью модели трения, изображенной на рис. 19. Объединив для удобства данную модель в один блок, промоделируем ее с помощью схемы, представленной на рис. 21. Результаты моделирования представлены на рис. 22.

Рис. 21. Структурная схема для моделирования смещения перед скольжением.

Полученная характеристика F(x) полностью совпадает с графиком, представленным в [9]. Отсюда следует, что даже тогда, когда приложенная сила не превышает уровень статического трения, есть микроскопическое смещение. Из графика зависимости F() видно, что при реверсе скорости, реверс силы трения имеет некоторую задержку. Данное свойство называют эффектом Даля (Dahl, 1968). Если на вход модели подать сигнал Fa > F, то тело выйдет из состояния малых смещений и произойдет резкое увеличение скорости.

Теперь покажем, как моделируется фрикционное запаздывание.

Сущность фрикционного запаздывания заключается в том, что трение и скорость связаны гистерезисной зависимостью. Сила трения меньше для убывающей скорости и больше для возрастающей скорости. Ширина петли гистерезиса растет с увеличением частоты изменения скорости. На рис. 23, б) представлены результаты моделирования с помощью модели, изображенной на рис. 23, а). На вход модели подана синусоида амплитудой Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

А=0.0007 м/с с постоянной составляющей, равной 0.001 м/с. Это означает, что реверса механизма нет. Моделирование выполнено для трех частот: = 1 рад/с, 2 = 10 рад/с и 3 = 25 рад/с, соответствующих опыту, проведенному в [5].

Рис. 22. Результаты моделирования смещения перед скольжением.

Рис. 23. Структурная схема и результаты моделирования Проведем опыт по определению влияния фрикционного запаздывания на характеристику трения при прохождении скорости через нулевое значение, т. е. при реверсе. Подадим на вход модели трения синусоиду с нулевой постоянной составляющей, амплитудой А=0.08 м/с, и частотой 1 = 10 рад/с и 2 = 20 рад/с. Результаты моделирования представлены на рисунке 24.

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

На основании эксперимента с движением «стик-слип» можно исследовать заключительное динамическое свойство данной модели, представленное в [9]: переменная сила трения при отрыве. Данное свойство характеризуется тем, что максимальная сила трения при переходе тела из состояния покоя (микроскопического движения) в состояние скольжения уменьшается с увеличением скорости нарастания силы, приложенной к телу. Исследуем данное свойство с помощью модели, представленной на рис. 25, а). Выполним моделирование с целью определения переменной силы трения при отрыве для четырех значений скорости нарастания приложенной силы Fa: 1 Н/c, 5 Н/c, 15 Н/c и 50 Н/c. Графики, полученные в результате моделирования, изображены на рис. 25, б) Из полученных графиков видно, что при увеличении скорости нарастания силы, приложенной к телу, сила трения при отрыве (т. е. статическое трение) снижается. Пределом для снижения является уровень Кулонова трения; в данном случае FC = 1 Н.

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

Рис. 25. Структурная схема и результаты эксперимента по определению переменной силы трения при отрыве.

3.2. Эласто-пластическая динамическая модель трения В последние годы были созданы еще две динамические модели трения, являющиеся модификациями модели ЛуГре и призванные улучшить ее работу в режиме предварительного смещения. Это эласто-пластическая модель и модель «Лойвен» (по названию университета, где работают ее авторы — Leuven).

Эласто-пластическая модель [16] устраняет недостаток, который существует как у статических моделей трения, так и у модели ЛуГре. Этим недостатком является так называемое «сползание». Эффект «сползания»

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

Для исследования подверженности моделей трения сползанию может быть проведен следующий тест, результаты которого представлены на рис. 26. К телу, покоящемуся на поверхности, была приложена сила Fa, которая сначала превысила силу трения покоя FS, затем была уменьшена и подавалась в виде синусоиды с постоянной составляющей, и в конце была снята (рис. 26, а). Для моделирования были использованы статическая модель с плавным переходом (рис. 26, в), модель ЛуГре (рис. 26, г) и эластопластическая модель (рис. 26, б). Очевидно, что и статическая модель, и модель ЛуГре «ползут», в отличие от эласто-пластической модели, совершающей лишь микроколебания, соответствующие упругому контакту.

В некоторых статьях показано, что модель ЛуГре, кроме того, неточно отражает гистерезисную зависимость при предварительном смещении.

Модель Лойвен [17] решает эту проблему, но она имеет достаточно сложное описание, предусматривающее использование большого количества внутренних состояний, тогда как модель ЛуГре и эласто-пластическая моСекция 5. Моделирование в Simulink дель имеют лишь одну переменную состояния. Таким образом, на сегодняшний день эласто-пластическая модель является наиболее полной моделью трения, которая по сложности сравнима с моделью ЛуГре.

4. Практическое применение библиотеки моделей трения Библиотека использовалась при идентификации параметров динамической модели трения [18] и разработке алгоритмов настройки компенсаторов его негативного влияния на станочный электропривод [12], [19]. Для проведения экспериментов использовался электропривод с винтовой передачей (ISW), оснащенный датчиками положения и ускорения (рис. 27). В ISW, руководимом профессором G. Pritschow, продолжаются активные исследования станочного привода, в частности, в области трения [20], [21].

Сбор данных и обработка измерений осуществлялся с помощью интерфейса фирмы dSPACE в совокупности со средствами RTW пакета MATLAB/Simulink при шаге дискретизации Td = 125 µs.

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

Рис. 27. Экспериментальный станочный электропривод (ISW).

Частотная характеристика механической части привода была экспериментально определена и рассчитана на основе дискретного преобразования Фурье, алгоритм которого реализован в МАТЛАБ в виде функции fft.

Вид ЛАЧХ представлен на рис. 28. Как видно из рис. 28, частотная характеристика механики привода имеет две доминирующие резонансные частоты (f1 = 117 Hz, коэффициент демпфирования D1 = 0.055, f2 = 349 Hz, коэффициент демпфирования D2 = 0.05) и ее аппроксимация может быть представлена передаточной функцией четвертого порядка. Экспериментально полученные зависимости момента трения от скорости представлены на рис. 29.

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

На рис. 30 изображены подобные зависимости, полученные в результате моделирования с использованием модели и следующих оценок параметров: F C + = F C = 1.6 Нм, F S + = F S = 1.85 Нм, S + = 0.004 м/с, S = 0.005 м/с, 2+ = 17 Нс/м, 2 = 8 Нс/м.

Экспериментально определенная кривая в режиме предварительного смещения (кривая Даля — зависимость момента трения от положения для малых смещений) представлена на рис. 31.

Рис. 28. ЛАЧХ механической части привода.

Рис. 29. Характеристика трения (эксперимент).

Рис. 30. Характеристика трения (моделирование).

Рис. 31. Характеристика трения в режиме предварительного смещения Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

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

Литература 1. Armstrong-Hlouvry B., Dupont P., Canudas De Wit C. A survey of models, analysis tools and compensation methods for the control of machines with friction // Automatica.— 1994.— V.30.— N.7.— P.1083–1138.

2. Первозванский А. А. Системы с разрывными нелинейностями при высокочастотных возмущениях // Автоматика и телемеханика.— 2000.— №7.— С.44–54.

3. Огурцов А. И. Модель плоского возмущенного движения ползуна с учетом нелинейности подъемной силы // СТИН.— 2000.— №7.— С.11–13.

4. Hess D. P., Soom A. Friction at a lubricated line contact operating at oscillating sliding velocities // J. Tribology.— 1990.— V.112.— P.147–152.

5. Андрющенко В. А. Следящие системы автоматизированного сборочного оборудования.— Л.: Машиностроение. Ленингр. отд-ние, 1979.— 246 с.

6. Крагельский И. В., Гитис Н. В. Фрикционные автоколебания.— М. Наука, 1987.— 183 с.

7. Заковоротный В. Л. Нелинейная трибомеханика.— Ростов-на-Дону:

Издательский центр ДГТУ, 2000.— 293 с.

8. Kozlowski K. R., Dutkiewicz P. Experimental identification of robot and load dynamics // IFAC 13th Triennial World Congress.— San Francisco, 1996.— P.397–402.

9. Canudas De Wit С., Olsson H., strm K. J., Lischinsky P. A new model for control of systems with friction // IEEE Transactions on Automatic Control.— 1995.— V.40.— N.3.— P.419–424.

10. Техническая кибернетика / Под ред. В.В. Солодовникова, кн. 3.— М.:

Машиностроение, 1967.— 770 с.

11. Kern R. H., Gao C.-T., Nitsche R. Reibkraftkompensation mittels FuzzyLogik // Automatisierungstechnischepraxis.— 1995.— V.37.— P.50–60.

12. Хлебалин Н. А., Костиков А. Ю. Моделирование трения и компенсация его негативного влияния в процедуре автоматизированной настройки регулирующего устройства технологической машины // Тр. международ. конф. «Идентификация систем и задачи управления» SICPRO’2003.

— М.: ИПУ РАН, 2003.— С. 1821–1849 (CD ISBN 5-201-14948-0).

13. Armstrong B. Challendges to Systematically Engineered Friction Compensation // Proc. of IFAC Workshop on Motion Control. Munich.— 1995.— P.21–30.

14. Canudas-De-Wit С. Comments on «A New Model for Control of Systems with Friction // IEEE Transactions on Automatic Control.— 1998.— V.43.— N.8.— P.1189–1190.

15. Пожбелко В. И. Законы предельного трения // Вестник Российской академии транспорта (Уральское межрег. отд.).— 1999.— №2.— С.226– 16. Dupont P., Hayward V., Armstrong B., Altpeter F. Single state elasto-plastic friction models // IEEE Transactions on Automatic Control.— V.47.— N.5.— May 2002.— P.787–792.

17. Swevers J., Al-Bender F., Ganseman C.G., Prajogo T. An integrated friction model structure with improved presliding behavior for accurate friction compensation // IEEE Transactions on Automatic Control.— V.45.— N.4.— April 2000.— P.675–686.

18. Костиков А. Ю, Хлебалин Н.А. Идентификация параметров динамической модели трения в электроприводах технологических машин// Тр.

международ. конф. «Идентификация систем и задачи управления»

SICPRO’2004.— М.: ИПУ РАН, 2004.— С.571–588 (CD ISBN 5-201Хлебалин Н. А. Идентификация параметров механической системы станочного электропривода с целью расчета регуляторов и диагностики // Тр. международ. конф. «Идентификация систем и задачи управления»

SICPRO’2000.— М.:ИПУ РАН.—С.1696–1709(CD ISBN 5-201-09605-0).

20. Pritschow G., Ketterer G. Rechnergestuetzte Identifikation von Reibkennlinien an elastisch gekoppelten Bewegungsachsen // Antriebstechnik.— 1993.— V.32.— N.9.— P.67–72.

21. Augsten G., Schmid D. Einfluss von Spiel und Reibung auf die Konturfehler bahngesteuerter Werkzeugmaschinen // Steuerungstechnik.— 1969.— V.2.— N.3.— P.103–108.

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

УДК

МОДЕЛИРОВАНИЕ ЭКОНОМИКИ В MATLAB/SIMULINK

Государственный университет управления, Москва, Excel прекрасный, почти святой, инструмент экономиста. Но для анализа, прогноза и планирования экономической динамики он слабо приспособлен. А ведь вся экономическая теория рассматривает в основном именно динамические проблемы, гипотезы и закономерности.

В то же время в математических и технических вузах РФ широко преподается математический программный инструмент MATLAB. Для него как надстройки (Toolboxes) разработаны многие спецприложения для анализа технических систем управления. Он также предоставляет экономистам финансовый пакет FinancialToolbox, связь с Excel — ExcelLink, связь с Word — Notebook.

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

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

Далее представлены три простейшие модели экономики в Simulink.

Модели преподаются в нескольких ВУЗах и позволяют сделать содержательные выводы.

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

Проблема: теория и практика не знают величину приемлимой для всех налоговой ставки. Задача: обосновать величину налоговой ставки.

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

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

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

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

Государство объявляет ставку налога на прибыль и получает от фирм средства в бюджет. Фирмы обладают собственным капиталом, производят прибыль, отчисляют по налоговой ставке средства в бюджет. Постналоговая прибыль как нераспределенная прибыль полностью включается в собственный капитал фирмы. Дивиденды не выплачиваются, никаких других отчислений от прибыли не производится. Вся прибыль распределяется только на два потока: в бюджет, а остаток в собственный капитал банка.

Математическая модель.

Сумма налоговых поступлений в бюджет за моделируемый период представлена формулой:

где BDt — сумма, поступивших в бюджет средств от начала моделирования к концу года t, руб.; PRFt — доналоговая прибыль (profit), полученная предприятием за год t, руб/год.; TXRT — ставка налога на прибыль (tax rate); t — время, год. Для запасов это момент — конец года, для потоков это интервал времени определенного года; tb — начальный (begin) интервал моделирования; tf — последний (final) год моделирования.

Капитализируемый предприятием за период моделировния остаток прибыли Прибыль за t год где RN — рентабельность капитала предприятия. Задается как параметр предприятия, исходное данное.

Компьютерная модель в программе Simulink. Компьютерная Simulink-модель представляется в виде блок-схемы, содержащей типовые функциональные блоки систем управления и управляемых объектов. В блоки включены компьютерные программы, вычисляющие математические функции. Значки на блоках представляют формулы аналитических Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

выражений передаточных функций как отношение выходной информации к входной. Окно модели представлено на рис. 1.1.

Рис. 1.1. Модель для определения оптимальной ставки налогообложения На схеме левый блок с именем «Бизнес» представляет накопитель собственного капитала предприятия, из библиотеки элеметов он взят как блок «Интегратор дискретного времени». На вход блока поступает поток капиталовложений CapF. Это постналоговая, нераспределенная прибыль.

Она аккумулируется, бизнесом и увеличивает его собственный капитал.

Выход блока — это величина капитала бизнеса CapS.

Справа от бизнеса блок умножения (с крестиком) производит поток прибыли PrfF как произведение капитала предприятия на рентабельность Rntb.

Рентабельность задается библиотечным блоком «Константа» с именем Rntb. Под ним как комментарий задан вектор плана экспериментов по фактору рентабельности [0.2:0.2:0:8]. В матричной лаборатории MATLAB это означает, что мы начнем имитационные эксперименты для рентабельности равной 20% с шагом в 20% и последним экспериментом для 80%.

Для автоматизации экспериментов константу рентабельности мы можем заменить на переменную и управлять ее значениями из программного файла MATLAB.

Ниже блока рентабельности расположен блок Scope1 — это графопостроитель для отображения переменной капитала предприятия CapS.

Следующий блок умножения (с крестиком) выдает произведение потока прибыли на налоговую ставку TaxRate. Это поток отчислений от прибыли в госбюджет TxF.

Ставка налога задается библиотечным блоком «Константа» с именем TaxRate. Под ним как комментарий задан вектор плана экспериментов по фактору налоговая ставка [0:0.1:1:0]. Это означает, что мы начнем имитационные эксперименты для ставки равной 0% с шагом в 10% и последним экспериментом для 100%. Для автоматизации экспериментов константу налоговой ставки мы можем заменить на переменную и управлять ее значениями из программного файла MATLAB.

Блок Госбюджет представлен интегратором. Он аккумулирует налоговые поступления TxF за период моделирования в виде переменной BdjS.

Справа от него блок Scope строит график накопления средств от налога в бюджете. Блок Display отображает числовые значения BdjS.

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

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

Средства управления экспериментом. Средства отображения информации о показателях экспериметов — это графопостроители Scope и индикаторы чисел Display.

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

Программа управления экспериментами %Optimal profit tax rate simulation %File: C:\Csr_MtLb\TxRt\TaxRate_DscM.m and TaxRate_Dsc.mdl %================================================= path(path,'C:\Csr_MtLb\TxRt') %Set model Path %================================================= open_system('TaxRate_Dsc') %Load TaxRate_Dsc.mdl TaxRate=[0:0.05:0.7] %План-вектор эксперимента по ставке налога for Rntb = 0.2:0.2:1 %Цикл и план-вектор по рентабельности sim('TaxRate_Dsc') %Run model plot(TaxRate, ScopeData(end,2:end)) %Чертить график поступлений в бюджет hold on %Разрешить дополнение графика кривыми grid %Чертить сетку end hold off %Запретить дополнение графика %======= end TaxRate_DscM.m program =============== В m-файле программы за знаком процента всегда идут поясняющие комментарии. Они не являются командами и компьютером не исполняются.

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

В первой строке программы дается ее назначение или смысловое название. Во второй строке полное имя.m-файла содержащего нашу программу для управления экспериментами над Simulink моделью и имя файла Simulink модели с расширением.mdl.

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

Шестая строка командой open_system загружает с диска модель в оперативную память.

Седьмая строка присваивает переменной модели TaxRate вектор плана экспериментов по налоговой ставке.

В строках с 8 по 14 выполняется for цикл для проведения экспериментов при различных величинах рентабельности бизнеса.

В девятой строке командой sim запускается модель и начинается моделирование, имитация налогового взаимодействия государства и предприятия.

После окончания имитации команда plot чертит один график рис.

1.3, используя данные рабочего (work space) пространства MATLAB, записанные туда графопостроителем Scope. Оператор hold on разрешает дополнять рисунок графиками кривых, рассчитанными для следующих значений циклов рентабельности.

Порядок выполнения работы. Вид экрана лабораторной модели представлен на рис. 1.2.

Рис. 1.2. Вид экрана лабораторной модели для исследования оптимальной ставки налогообложения.

Слева расположено окно редактора программ с частичной видимостью m-файла программы управления двухфакторным экспериментом.

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

Задание №1. Однофакторный имитационный эксперимент. Исследовать зависимость налоговых поступлений в бюджет за конкретный период времени от величины налоговой ставки на прибыль предприятий.

Запустив модель из меню Simulation, наблюдаем в окнах Scope изменение показателей предприятий и бюджета во времени: рост поступлений прибыли, отчислений по налогу в бюджет и капитализацию нераспределенной прибыли бизнесом. Устанавливая различные ставки налога каждый раз прогоняем модель.

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

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

В командном окне MATLAB откроем вышеописанный файл TaxRate_DscM.m. Файл откроется в окне редактора. Этот файл будет управлять двухфакторным экспериментом. Запустим программу командой меню Tools > Run. В результате моделирования MATLAB построит графики зависимости поступлений в бюджет от налоговой ставки и рентабельности предприятий, рис. 1.3.

Анализ результатов. Чем выше рентабельность предприятия, тем ярче выражена оптимальная ставка налогообложения. С ростом рентабельности оптимальная ставка уменьшается (сдвигается влево), стремясь к фиксированной величине, на наших графиках, примерно, к 23%.

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

Рис. 1.3. Зависимость поступлений в бюджет от налоговой ставки для предприятий различной рентабельности.

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

2. Равновесие на конкурентном рынке Определение проблемы.Основоположник ценовой теории Альфред Маршалл (1842–1924) полагал, что большинство экономических процессов можно объяснить в терминах равновесной рыночной цены. Цена устанавливается при взаимодействии спроса и предложения. Теории хороши во всех учебниках, но и через 80 лет после Альфреда богатейшая и грамотнейшая OPEC, регулируя предложение, не может установить приемлимые потребителям и производителям цены на нефть.

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

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

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

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

Функция зависимости спроса от цены где Dmd — спрос(demand) за текущий интервал времени; D0 — спрос при нулевой цене; Kd — крутизна линии спроса; Prc — цена (price) товара.

Линия зависимости предложения от цены где Spl — предложение(supply) за текущий интервал времени; S0 — предложение при нулевой цене; Ks — крутизна линии спроса; Prc — цена (price) товара.

Имитационная модель в Simulink. Блок-схема имитационной модели представлена на рис. 2.1.

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

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

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

Обозначения и параметры блока на схеме следующие:

Блок Лаг имитирует запаздывание поставщика на рынке. Продавец поставляет товар в количестве Spl, определенном на основе цен прошлого интервала времени.

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

Параметры блока одинаковы с параметрами блока Спрос DmdFn.

Средства управления экспериментом.

Внутримодельные средства.

Кроме вышеописанной экономической части модели в ней присутствуют блоки управления экспериментом. Отображение результатов моделирования осуществляется блоками Scope. На схеме модели они расположены справа. Блоки строят графики изменения во времени спроса, предложения и цены. Слева блоки константы, интегратора и усилителя задают значения цены для построения графиков статических характеристик функций спроса и предложения.

Ключи предназначены для переключения режимов моделирования.

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

Внешние средства — программа MATLAB. Автоматизация управления экспериментом осуществляется с помощью программы MATLAB.

%Market equlibrim price simulation %File: C:\Csr_MtLb\DmdSplEquM.m %1.Simulate static functions, plot its %2.Simulate price dynamic, plot price Web graphics %================================================= %Set model Path path(path,'C:\Csr_MtLb\MrktEqlPrc') %================================================= %1.Simulate static functions, plot its %Load and Run DmdSplEqu.mdl open_system('DmdSplEqu') sim('DmdSplEqu')%Write Vars into WS from Scope %Plot Static features plot(ScopeData(:,2),ScopeData(:,3:4)) hold on grid pause(5)%Пауза для переключения в режим имитации %двойным щелчком мыши над ключами %================================================= %2.Simulate price dynamic sim('DmdSplEqu') %================================================= %3. plot price Web graphics for i = 2:11 %Цикл черчения паутины движения к равновесию line([ScopeData(i-1,2) ScopeData(i,2)],[ScopeData(i,4) ScopeData(i,4)]) line([ScopeData(i,2) ScopeData(i,2)],[ScopeData(i,4) ScopeData(i+1,4)]) end hold off %================================================= Описание программы. Строки со знаком % представляют комментарий к программе и не выполняются компьютером. Командой path устанавливается путь к модели в файловой системе. Оператор open system загружает блок-модель в Simulink.

Функция sim запускает модель. Plot чертит крест графиков спроса и предложения. Hold on разрешает дополнять рисунок новыми графиками.

Grid рисует для графиков масштабную сетку. На этом первый этап заканчивается. Функция pause останавливает процесс моделирования на 5 секунд для установки ключей двойным щелчком мыши в нижнее положение.

Начинается второй этап — моделирование переходного процесса к рыночному равновесию. Sim повторно запускает модель. Идет имитация.

Scopes в своих окнах чертят графики показателей и пишут их значения в рабочее пространство MATLAB workspace.

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

Общий вид окон эксперимента представлен на рис. 2.2.

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

Рис. 2.2. Вид экрана при выполнении лабораторных работ по анализу процессов рыночного равновесия.

Слева вверху расположено окно редактора/отладчика MATLAB c мфайлом управления экспериментом. Справа вверху окно модели Simulink.

Внизу слева окно Scope, справа Scope1 с графиками изменения цены, спроса и предложения во времени. В центре график паутины переходного процесса к равновесию на конкурентном рынке, построенный программой MATLAB.

Задание №1. Изучить переходный процесс к рыночному равновесию. Из командного окна MATLAB студент открывает m-файл DmdSplEquM.m. Файл отображается в окне редактора. Запускаем файл из меню Tools>Run. Программа загружает и стартует модель. В окне Simulink переводим ключи в нижнее положение. Наблюдаем графики изменения показателей во времени в окнах Scope и паутину движения цены в окне Figure. Графики Scope1 представлены на рис. 2.3.

Графики Scope представлены на рис. 2.4.

Рис. 2.3. Изменение цены, спроса-предложения во времени.

Рис. 2.4. Изменение цены, спроса и предложения в окне Scope.

График паутинообразной модели движения к рыночному равновесию представлен на рис. 2.5.

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

Рис. 2.5. График паутины движения цены к равновесию на рынке.

Задание №2. Изучить влияние смещения линий спроса и предложения на рыночное равновесие. Для смещения спроса студент меняет значения параметра D0 в блоках DmdFn и SplFn1 и наблюдает движение цены к новому равновесию.

Для смещения предложения студент меняет значения параметра S0 в блоке SplFn и наблюдает движение к новому равновесию.

Задание №3. Изучить влияние крутизны линий спроса и предложения на рыночное равновесие. Для изменения крутизны спроса студент меняет значения параметра Kd в блоках DmdFn и SplFn1 и наблюдает движение цены к новому равновесию.

Для изменения крутизны предложения студент меняет значения параметра Ks в блоке SplFn и наблюдает движение к новому равновесию.

При увеличении Ks колебания показателей в системе возрастают.

При Ks = Kd параметры колебаний не изменяются во времени. Устойчивое равновесие не достигается. На экране паутина превращается в прямоугольную рамку, т. е. система находится в режиме устойчивых колебаний.

При Ks > Kd система неустойчива, процесс расходится, паутина расширяется и уходит из окна. Это противоречит реальной экономике и подтверждает, что модель груба.

3.Циклы и кризисы Определение проблемы. Мировую экономику и экономику отдельных стран периодически потрясают кризисы или циклы роста и падений.

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

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

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

Например, в 1962 году выпускали шесть миллионов швейных машин, а через год только один миллион. Шесть заводов уничтожили из-за перепроизводства. То же было с комбайнами, литейным оборудованием и др. Даже по очковым линзам, где мода и технология не менялись 300 лет, каждые 10 лет создавались и уничтожались заводы и цеха. На этих заводах менялся только привод станков для шлифовки линз: ручной, ножной, конский, водяной, паровой, электрический.

Причина циклов и кризисов остается во многом не ясной. В этом проблема.

Цель работы. Исследовать причинно-следственный механизм возникновения циклов и кризисов перепроизводства.

Лабораторная модель. Имея Simulink с типовыми библиотечными блоками мы можем не выписывать математическую модель — каждая элементарная модель уже имеет программу, привязанную к блоку. Блоксхемная имитационная модель Simulink представлена на рис. 3.1.

Рис. 3.1. Модель для анализа циклов перепроизводства.

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

Последние 70 лет кризисы перепроизводства, как правило, начинались в автомобильной промышленности и через 5-10 лет после окончания войн. Поэтому предположим, что промышленность выпускает оборудование в виде автомобилей.

Модель потребности в парке оборудования на схеме задается слева в верху блоком констант и интегратором. На выходе блоков линейная во времени переменная потребности населения в парке оборудования (автомобилей) с именем NeedPrk (необходимый парк). Ниже потребностей расположены четыре блока, отражающие движение парка оборудования: поступление оборудования в эксплуатацию от производителей (переменная Prdct), накопление их в парке интегратором c выходной переменной InPrk ( в парке). Блоки срок службы и выбытие задают переменную OutPrk (поток выбытия оборудования по ветхости, износу или моральному старению).

Нижний круглый блок вычитает из поступившего в парк оборудование выбывшее, создавая переменную RealPrk (реальное количество товаров, находящихся в эксплуатации). Верхний круглый блок сумматора вычитает из необходимого парка наличный парк создавая переменную текущего спроса СurDmd. Блоком Saturation (ограничитель) она обрезается снизу, реализуя традиционную для экономических задач неотрицательность переменных:

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

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

В данной модели смена параметров блоков производится вручную.

Диалоговое окно настройки параметров вызывается двойным щелчком мыши на блоке.

Полный набор окон эксперимента представлен на рис. 3.3.

Секция 5. Моделирование в Simulink Рис. 3.2. Simulink модель анализа циклов и кризисов с блоками графопостроителей Scope.

Рис. 3.3. Экран с полным набором окон эксперимента для исследования циклов и кризисов.

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

Слева вверху расположено окно библиотеки стандартных блоков Simulink. Справа окно модели. Слева внизу окно графопостроителя Scope1.

В центре командное окно MATLAB. Справа окно Scope.

Задание №1. Исследовать влияние производственного лага на устойчивость экономики. Студент на рабочем столе компьютера двойным щелчком мыши на иконке MATLAB загружает его. В командном окне вызывает Simulink и открывает файл ParkPrdct.mdl.

Запускает модель из меню Simulation>start. В окнах Scope наблюдает графики изменения показателей экономического развития, представленные на рис. 3.4 и 3.5.

Рис. 3.4. Графики экономической динамики производства, спроса и потребления, Рис. 3.5. Графики экономической динамики производства, спроса и потребления, Scope.

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

Задание №2. Исследовать влияние срока службы изделий на динамику производства. Студент произвольно меняет значения параметра задержки в блоке Срок службы, запускает модель, наблюдает графики показателей. С увеличением задержки возрастает амплитуда и период колебаний показателей экономической системы, т. е. снижается ее устойчивость, возрастает неустойчивость, возможность кризисов. Это характерно для отраслей, производящих товары долгосрочного пользования или длительного хранения. Для товаров разового потребления, не подлежащих длительному хранению, при нулевых задержках производства колебания не возникают. Предложение всегда равно спросу. Система устойчива.

Задание №3. Исследовать влияние начального дефицита на устойчивость производства. Под начальным дефицитом понимается разница между необходимым и реальным парком оборудования на момент начала моделирования.

Студент произвольно меняет величину начальных условий на интеграторе блока Поступление, изменяя величину начального реального парТруды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

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

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

Литература 1. Гайгер Л. Т. Макроэкономическая теория и переходная экономика.— М.: ИНФРА-М,1996.

2. Пиндайк Р., Рубинфельд Д. Микроэкономика.— М.: Дело, 1992.

3. Цисарь И. Ф Крыкин М. А. MATLAB_Simulink — лаборатория экономиста.— М.: Анкил, 2001.

4. Цисарь И. Ф. Лабораторные работы на персональном компьютере.— М.: Экзамен, 2002.

5. Цисарь И. Ф Нейман В. Г. Компьютерное моделирование экономики.— М.: Диалог-МИФИ, 2002.

6. Черник Д. Г. и др. Налоги.— М.: Финансы и статистика, 1997.

УДК 613.

МЕТОДИКА МОДЕЛИРОВАНИЯ УСТАНОВОК

ИНДУКЦИОННОГО НАГРЕВА В SIMULINK

Уральский государственный технический университет, Екатеринбург, e-mail:[email protected], [email protected] Российская электротехнологическая компания (РЭЛТЕК), Екатеринбург, В настоящее время существует несколько пакетов схемотехнического моделирования: OrCad, Micro Cap, Electronic Work Bench и т. п. Данные пакеты в основном ориентированы на разработку электронных схем, включая проектирование печатных плат, и не позволяют в полной мере проводить моделирование силовых электронных схем, особенно если необходимо выполнить расчеты преобразователя энергии совместно с электротехнической установкой или электрической машиной. Моделирование системы управления силовой части преобразователя в таких пакетах требует больших трудовых затрат, поскольку систему управления необходимо создавать из моделей реальных электронных компонентов, в то время как достаточно часто требуется лишь проверить алгоритм ее работы. И, наконец, модели полупроводниковых приборов, используемые в таких пакетах, являются очень подробными (моделируется плавная вольт-амперная характеристика, учитываются динамические свойства и т. д.). В большинстве случаев такая детализация не нужна, достаточно учесть лишь основные свойства полупроводниковых приборов. Например, для диода достаточно учесть свойство односторонней проводимости и внутреннее сопротивление. Излишне подробное описание полупроводниковых приборов приводит также к тому, что схема состоящая из нескольких приборов рассчитывается слишком долго.

В этом смысле инструмент моделирования динамических систем Simulink пакета MATLAB совместно с библиотекой электротехнических блоков Sim Power System (SPS) выгодно отличается от пакетов схемотехнического моделирования. Библиотека SPS содержит широкий набор моделей элементов электротехнических цепей, электрических машин и компонентов силовой электроники. При этом точность моделей такова, что основные свойства перечисленных объектов учитываются в достаточной мере. При моделировании электротехнологических установок силовая часть преобразователя энергии создается при помощи блоков библиотеки Sim Power System, а модель обрабатываемого объекта (тепловая часть модели) и система управления из блоков основной библиотеки Simulink.

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

Однако совместный расчет преобразователя и индуктора нагревательной установки во временной области требует больших затрат времени. Это связано с тем, что постоянные времени в тепловой части модели (нагреваемое тело) существенно больше, чем в электрической ее части (преобразователь энергии и система управления). Это означает, что тепловая часть модели должна рассчитываться с избыточно мелким шагом, а значит, время расчета всей установки будет существенно увеличено, поскольку общее время расчета определяется именно временем расчета тепловой части модели. Например, при частоте выходного напряжения преобразователя 1000 Гц и числе расчетных точек на один период напряжения равном 50, требуемый расчетный шаг электрической части модели будет составлять 210-5 с. В то же время, тепловая постоянная времени нагреваемого тела может составлять 10–100 с и, соответственно, шаг расчета тепловой части модели может быть значительно больше: 0.01–0.001 c. Данное соотношение между требуемыми расчетными шагами в тепловой и электрической частях модели означает, что на каждый расчетный шаг тепловой части модели должно приходиться 50-500 шагов в электрической части модели, что в такое же число раз увеличивает время расчета установки в целом.

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

1. Начало расчета.

2. Расчет электрической части модели с малым шагом на некотором интервале времени. Тепловая часть модели в это время не рассчитывается.

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

4. Передача значения выходной активной мощности преобразователя в тепловую часть модели с одновременным прекращением расчета электрической части модели и началом расчета тепловой части модели.

Расчетный шаг при этом существенно увеличивается.

5. Определение температуры нагреваемого тела в конце интервала расчета тепловой части модели для передачи в электрическую часть модели.

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

7. Уменьшение шага по времени до величины требуемой для расчета электрической части модели.

8. Переход к п. 2 алгоритма.

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

Общая схема модели показана на рис. 1.

Электрическая часть модели представлена подсистемой Electro, а тепловая подсистемой Teplo. Блок Pulse Generator задает прямоугольный логический сигнал, скважность которого определяет длительность расчета каждой из стадий расчета. При этом нулевое значение сигнала переключает схему на расчет электрической части модели (подсистема Electro), а единичное значение сигнала переводит схему на расчет тепловой части модели (подсистема Teplo). Подсистема Step Control Subsystem обеспечивает формирование сигналов управления запускающих на расчет электрическую и тепловую части модели. Схема данной подсистемы представлена на рис. 2.

Рис.2. Схема системы управления расчетом.

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

Входной прямоугольный сигнал через блок умножения Product проходит на первый выход подсистемы и обеспечивает включение тепловой части модели. Этот же, но инвертированный сигнал проходит на второй выход и обеспечивает, в свою очередь, запуск электрической части модели. Таким образом, электрическая и тепловая части модели работают в противофазе, а длительность работы каждой из частей модели определяется скважностью входного сигнала. Блок Timer задает увеличенный первый интервал расчета электрической части модели, для того, чтобы преобразователь энергии вышел на установившийся режим работы. Переключение шага расчета выполняют триггерные (срабатывающие по фронту сигнала) подсистемы Step Switch. Подсистема Step Switch (рис. 3) использует блок Constant для задания величины шага, а блок S-function выполняет установку величины шага в процессе расчета модели. На рис. 3 показано значение шага для режима расчета тепловой части модели.

Величина требуемого шага устанавливается блоком S-function с помощью команды вида:

set_param(имя_модели,'maxstep',num2str(u) ), где set_param — имя команды устанавливающей текущее значение максимального шага расчета, имя_модели — имя файла Simulink-модели, maxstep — имя устанавливаемого параметра (максимальный шаг расчета), u —числовое значение устанавливаемого шага.

Текст S-функции приведен ниже:

function [sys,x0,str,ts]=Set_Max_Step_Size(t,x,u,flag,m_name) switch flag, case 0, [sys,x0,str,ts]=mdlInitializeSizes(m_name);

case 3, sys=mdlOutputs(t,x,u,m_name);

case {1,2,4,9} sys=[];

otherwise error(['Unhandled flag = ',num2str(flag)]);

end function [sys,x0,str,ts]=mdlInitializeSizes(m_name) sizes = simsizes;

sizes.NumContStates = 0;

sizes.NumDiscStates = 0;

sizes.NumOutputs = 0;

sizes.NumInputs = 1;

sizes.DirFeedthrough = 1;

sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

x0 = [];

str = [];

ts = [-1];

function sys=mdlOutputs(t,x,u,m_name) set_param( m_name,'maxstep',num2str(u) );

sys=[];

В теле S-функции переменная m_name является именем рассчитываемой модели. Подробно методика разработки S-функций изложена в [1].

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

Схема электрической части модели показана на рис. 4.

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

Сигнал задания на величину сопротивления формируется из двух составляющих: первая — начальное значение сопротивления (при не нагретом рабочем теле индукционной установки), вторая — приращение сопротивления вследствие нагрева рабочего тела. Температура нагрева рабочего тела поступает на вход T_In из тепловой части модели. Активная мощность в индукционной установке формируется путем перемножения сигналов снимаемых с датчика тока и датчика напряжения. Апериодическое звено выполняет усреднение значения активной мощности и, затем сигнал пропорциональный мощности поступает в тепловую часть модели.

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

Рис. 5. Модель управляемого активного сопротивления.

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

На рис. 7 представлены графики сигналов моделируемой индукционной установки.

В примере на рис. 7, для повышения наглядности, время работы электрической части модели выбрано всего в три раза меньше, чем время работы тепловой части модели. В практических расчетах это соотношение может быть значительно увеличено. На графиках хорошо видно, что при Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

работе электрической части модели расчет ее тепловой части не производится и наоборот. Таким образом скорость расчета модели в целом оказывается существенно увеличенной.

Литература 1. Черных И. В. Simulink. Среда создания инженерных приложений.— M.: Диалог-МИФИ, 2004.— 496 с.

УДК 621.878.

МОДЕЛИРОВАНИЕ И ИССЛЕДОВАНИЕ РАБОЧЕГО

ПРОЦЕССА РЫХЛИТЕЛЬНОГО АГРЕГАТА С СИСТЕМОЙ

УПРАВЛЕНИЯ В СРЕДЕ SIMULINK

Сибирская государственная автомобильно-дорожная академия Омск, Рабочий процесс (РП) рыхлительного агрегата (РА) статического действия является сложной динамической системой. Первичный источник энергии, которым является двигатель внутреннего сгорания, работает в условиях постоянно изменяющейся нагрузки, носящей характер случайного процесса. Изменения нагрузки обусловлены неоднородностью разрабатываемого грунта, сколами грунта в процессе его разработки и случайными воздействиями от микрорельефа местности, по которой движется РА [6].

Использование системного подхода при анализе РП РА [2] позволило выявить следующие составляющие рабочего процесса (рис. 1): подсистема «микрорельеф», подсистема «остов — навесное оборудование», подсистема «рабочий орган — грунт», ходовое оборудование, гидропривод рабочего органа (РО), трансмиссия, двигатель внутреннего сгорания (ДВС) и система управления.

Рис. 1. Блок-схема рабочего процесса рыхлительного агрегата.

Составляющие РП РА взаимодействуют между собой через связи Pi.

Микрорельеф через подвеску воздействует на остов (P1), вызывая при этом неуправляемые перемещения РО (P2), а также на ходовое оборудование Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

(P11), вызывая дополнительные сопротивления движению (уклон местности, сопротивление перекатыванию). Ходовое оборудование, в свою очередь, воздействует на микрорельеф (P6), сглаживая его, и на разрабатываемый грунт (P5) (буксование). Сопротивление грунта воздействию РО, характеризуемое глубиной рыхления (Р3), вызывает сопротивление перемещению РА (Р4) и перемещения остова (Р7). Суммарный момент нагрузки (Р8) передается через трансмиссию (Р14) на ДВС (Р15), вызывая изменение угловой скорости его вала (Р10) и, соответственно, скорости движения РА, а следовательно и интенсивности взаимодействия ходового оборудования с микрорельефом и РО с разрабатываемым грунтом. Гидропривод РО изменяет глубину рыхления (Р11) и при работе вызывает дополнительное увеличение нагрузки на ДВС (Р16).

Система управления (СУ) получает информацию о нагрузке на двигатель (Р17) и о угловой скорости вала двигателя (Р10). На основе полученной информации СУ вырабатывает управляющие воздействия на ДВС (Р18), на трансмиссию (Р12) и на гидропривод рабочего органа (Р13), стремясь обеспечить оптимальную величину нагрузки на ДВС [7].

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

Микрорельеф описывается с помощью корреляционных функций, например микрорельеф целины [1]:

где — среднеквадратическое отклонение высотных отметок микрорельефа; А1, А2 — весовые коэффициенты составляющих; 1, 2 — коэффициенты затухания; — коэффициент периодичности.

Корреляционная функция (1) реализуется с помощью рекуррентного уравнения [1]:

y(n) = a01x(n) + b11y(n-1) + a02x(n) + a12x(n-1) + b12y(n-1) + b22y(n-2), (2) где а01, а02, а12, b11, b12, b22 — коэффициенты рекуррентного уравнения [1].

Сглаживающая способность гусениц учитывается уравнением:

где k = 0,5(Mc — 1); Mc — интервал усреднения; y(n) — ординаты несглаженного микрорельефа.

При описании колебаний остова, подвешенного на эластичной подвеске, использовался принцип Даламбера, при этом были получены следующие передаточные функции для вертикальных перемещений где y1, y2 — изменение вертикальной координаты под передней и задней полугусеницами; Fв — изменение вертикальной составляющей реакции грунта на РО; kiy, Tiy — соответственно коэффициенты передачи и постоянные времени, зависящие от коэффициентов жесткости и упругости элементов подвески.

Передаточные функции для угловых колебаний остова:

где МF — изменение момента, поворачивающего остов, от действия реакции грунта на РО.

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

где yро — изменение вертикальной координаты рабочего органа от неуправляемых перемещений остова; yвз — изменение вертикальной координаты под осью ведомой звездочки; kiро — коэффициенты усиления, зависящие от геометрических размеров РА; ро — время запаздывания.

где kб — коэффициент базы.

где Lб — длина базы.

Глубина рыхления складывается из заданной глубины рыхления hупр, которая формируется СУ, и неуправляемых перемещений РО yро.

Подсистема «РО — грунт» описывает зависимость сопротивления грунта рыхлению (горизонтальная Fг и вертикальная Fв составляющие) от глубины рыхления и скорости поступательного движения РА, при этом сопротивление грунта рыхлению складывается из тренда и флуктуаций, которые учитывают неоднородность грунта и его сколы в процессе рыхления [6].

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

Используя принцип Даламбера, поступательное движение РА можно описать уравнением [4]:

где Мв — момент на условном валу трактора; J3 — приведенный к валу трактора момент инерции РА от поступательного движения; Мf — момент сопротивления качению трактора; Мро — момент сопротивления, создаваемый РО; М — момент сопротивления, от уклона местности; 2 — угловая скорость вала трактора.

где Pf(vтр) — зависимость силы сопротивления качению трактора от скорости его поступательного движения; Мра — масса рыхлительного агрегата;

дв — уклон местности; б — коэффициент буксования, зависящий от силы тягового сопротивления; 1 — угловая скорость вала ДВС.

где rк — радиус ведущей звездочки; iтр — передаточное число трансмиссии.

Уравнение трансмиссии [6]:

где Мс — момент сопротивления, на валу ДВС; J2 — момент инерции вала трансмиссии; kтр — коэффициент, аппроксимирующий зависимость момента сопротивления сил трения в трансмиссии от угловой скорости вала трансмиссии.

Передаточная функция ДВС, учитывающая нелинейность механической характеристики [7]:

где Ср, Ск — коэффициенты передачи для регуляторной и корректорной ветвей механической характеристики, соответственно; Тр, Тк — постоянные времени для регуляторной и корректорной ветвей соответственно;

Мн — номинальный крутящий момент ДВС.

Передаточные функции гидропривода (ГП) [5]:

где Lшт — изменение длины штока гидроцилиндра; Мгп — дополнительный момент сопротивления, от включения гидропривода; U — сигнал управления гидроприводом; kгп1 — коэффициент передачи ГП, характеризующий скорость движения штока гидроцилиндра; kгп2 — коэффициент передачи ГП, характеризующий максимальное значение дополнительного момента сопротивления; Тгп — постоянная времени, характеризующая инерционность ГП; гп — время технического запаздывания ГП.

Обобщенная математическая модель рабочего процесса рыхлительного агрегата и СУ (рис. 2), составлена на основе зависимостей (2) (22), подсистема «Навесное оборудование» описывает зависимость перемещений РО от перемещений штока гидроцилиндра.

До настоящего времени в системах управления землеройными машинами в основном использовался способ управления, суть которого заключался в том, что информационными параметрами являлись выходные параметры рабочего процесса (вертикальная координата РО, угловая скорость вала ДВС и др.), которые после соответствующей обработки СУ формировали управляющее воздействие на гидропривод рабочего оборудования. Известны системы управления, построенные по принципу управления по возмущению, суть которых заключалась в измерении возмущающего воздействия или его отклонения. По величине этого возмущения формировалось управляющее воздействие (на двигатель или РО), и тем самым компенсировалась основная помеха, приводящая к потере оборотов двигателя, а следовательно и скорости машины [3, 7].

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

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

Информационными параметрами предложенной системы управления являются: угловая скорость вала ДВС (первый контур) и отклонение момента сопротивления на валу ДВС (второй контур), которые в соответствии с предложенным алгоритмом управления формируют управляющие воздействия, соответственно на рейку топливного насоса и на гидропривод РО, стабилизируя глубину рыхления и момент сопротивления на валу ДВС.

Труды II научной конференции «Проектирование инженерных и научных приложений в среде MATLAB»

Рис. 2. Структурная схема рабочего процесса рыхлительного агрегата с системой Управляющее воздействие на ГП РО формируется с помощью порогового элемента с регулируемыми зонами нечувствительности, статическая характеристика которого представлена на рис. 3. Максимальные отрицательное b1 и положительное b2 отклонения момента сопротивления на валу ДВС (левая и правая границы зоны нечувствительности) не равны (b1 b2), вследствие нелинейности механической характеристики ДВС. Зона нечувствительности c предназначена для нейтрализации дополнительного момента сопротивления, возникающего при включении ГП.

Рис. 3. Статическая характеристика порогового элемента СУ:

b1 — нижний порог срабатывания для корректорной ветви; b2 — верхний порог срабатывания для регуляторной ветви; с — порог срабатывания на отключение.

Функцией задатчика является задание оптимальной, для данных грунтовых условий и параметров микрорельефа, глубины рыхления hзад и значения крутящего момента на валу ДВС Мсзад, которое должна поддерживать СУ.

Представленная математическая модель РП РА с СУ была реализована в среде Simulink. Исследования РП РА без СУ позволили выявить зависимости значений оптимальной глубины рыхления, и соответствующего крутящего момента на валу ДВС, от характеристик разрабатываемого грунта и микрорельефа по которому движется РА.

Дальнейшие исследования направлены на выявления оптимальных параметров ГП РО и СУ, а также составления методики выбора этих параметров, в зависимости от характеристик разрабатываемого грунта и РА.

Литература 1. Глушец В. А. Математическая модель процесса взаимодействия гусеничного ходового оборудования землеройно-транспортных машин с разрабатываемым грунтом / Межвуз. сб. тр. ученых, аспирантов и студентов.— Омск: СибАДИ, 2004.— Вып. 1.— Ч.1.— С.152–158.

2. Глушец В. А. Математическая модель рабочего процесса рыхлительного агрегата // Дорожно-транспортный комплекс, экономика, экология, строительство и архитектура: Материалы Международ. науч.-практ.

конф. Кн. 2.— Омск: Изд-во СибАДИ, 2003.— С.250–253.

3. Ксеневич И. П., Тарасик В. П. Системы автоматического управления ступенчатыми трансмиссиями тракторов.— М.: Машиностроение, 1979.

4. Кутьков Г. М. Тяговая динамика тракторов.— М.: Машиностроение, 1980. — 215 с.

5. Титенко В. В. Повышение точности планировочных работ автогрейдерами с дополнительными опорными элементами рабочего органа:

Дис. … канд. техн. наук

.— Омск, СибАДИ, 1997.— 172 с.

6. Щербаков В. С., Глушец В. А. Математическая модель силового воздействия грунта на рабочий орган рыхлительного агрегата // Сб. науч. тр.

№4. Юбилейный.— Омск: Изд-во СибАДИ, 2003.— С.72–75.

7. Щербаков В. С., Глушец В. А. Система управления рыхлительным агрегатом для разработки мерзлых грунтов в условиях Севера // Проблемы создания и эксплуатации автомобилей, специальных и технологических машин в условиях Сибири и крайнего севера: Материалы Международ.

науч.-техн. конф.— Омск: Изд-во «ЛЕО», 2004.— С.219–220.