Численное моделирование процесса движения щита микротоннелепроходческого комплекса Текст научной статьи по специальности «Механика и машиностроение»

ISSN 0321-2653 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION. TECHNICAL SCIENCE. 2019. No 4

УДК 681.511 DOI: 10.17213/0321-2653-2019-4-44-50

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА ДВИЖЕНИЯ ЩИТА МИКРОТОННЕЛЕПРОХОДЧЕСКОГО КОМПЛЕКСА*

© 2019 г. А.В. Батюков1, Н.А. Глебов1, А.В. Павленко1, М.А. Земляной2, А.В. Живодерников1

1Южно-Российский государственный политехнический университет (НПИ) имени М.И. Платова, г. Новочеркасск, Россия,

2ООО ПСК «Гидрострой» г. Новороссийск, Россия

NUMERICAL SIMULATION OF THE MICROTUNNELING SHIELD MOVEMENT PROCESS

A.V. Batyukov1, N.A. Glebov1, A.V. Pavlenko1, M.A. Zemlyanoy2, A.V. Zhivodernikov1

1Platov South-Russian State Polytechnic University (NPI), Novocherkassk, Russia, 2LLS DCC «Gidrostroy», Novorossiysk, Russia

Батюков Александр Владимирович - инженер, НИИ Электромеханики, Южно-Российский государственный политехнический университет (НПИ) имени М.И. Платова, г. Новочеркасск, Россия. E-mail: Alexbatyukov@gmail.com

Глебов Николай Алексеевич — д-р техн. наук, профессор, кафедра «Мехатроника и гидропневмоавтоматика», ЮжноРоссийский государственный политехнический университет (НПИ) имени М.И. Платова, г. Новочеркасск, Россия. E-mail: aprim. srstu@mail.ru

Павленко Александр Валентинович — д-р техн. наук, профессор, зав. кафедрой «Электромеханика и электрические аппараты», Южно-Российский государственный политехнический университет (НПИ) имени М.И. Платова, г. Новочеркасск, Россия. E-mail: rn6lde@mail.ru

Земляной Михаил Александрович - канд. техн. наук, директор, ООО ПСК «Гидрострой», г. Новороссийск, Россия. E-mail: p218_gidrostroy@mail.ru

Живодерников Андрей Вячеславович — мл. науч. сотрудник, НИИ Электромеханики, Южно-Российский государственный политехнический университет (НПИ) имени М.И. Платова, г. Новочеркасск, Россия.

Batyukov Alexander Vladimirovich — Engineer, Engineer of the Research Institute of Electromechanics, Platov South-Russian State Polytechnic University (NPI), Novocherkassk, Russia. E-mail: Alexbatyukov@gmail.com

Glebov Nikolay Alexeevich - Doctor of Technical Sciences, Professor, Department «Mechatronics and Hydropneumoauto-matics», Platov South-Russian State Polytechnic University (NPI), Novocherkassk, Russia. E-mail: aprim.srstu@mail.ru

Pavlenko Alexander Valentinovich - Doctor of Technical Sciences, Professor, Head of Department «Electromecanics and Electric Devices», Platov South-Russian State Polytechnic University (NPI), Novocherkassk, Russia. E-mail: rn6lde@mail.ru

Zemlyanoy Mikhail Alexandrovich - Candidate of Technical Sciences, Director, LLS DCC «Gidrostroy», Novorossiysk, Russia. E-mail: p218_gidrostroy@mail.ru

Zhivodernikov Andrey Vyacheslavovich - Junior Researcher, Engineer of the Research Institute of Electromechanics, Platov South-Russian State Polytechnic University (NPI), Novocherkassk, Russia.

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

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

НИОКТР выполняется в ЮРГПУ(НПИ) при финансовой поддержке Министерства науки и высшего образования Российской Федерации в рамках Соглашения о предоставлении субсидии № 074-11-2018-010 от «5» июня 2018 г. по теме: "Создание высокотехнологичного производства импортозамещающего горно-проходческого оборудования, оснащенного интеллектуальными системами управления, для освоения подземного пространства городов".

ISSN 0321-2653 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION. TECHNICAL SCIENCE. 2019. No 4

To solve the problem of numerical simulation of the motion of the tunnel shield of the microtunneling complex, the article discusses the principles of construction and numerical simulation of an observer that allows to evaluate the coordinates of the tunnel shield inaccessible to measurement, and methods for creating the optimal controller in the automatic motion control system. The numerical model was built on the basis of the linearized equations of the tunneling shield's motion and the synthesized optimal controller. As a result of numerical mo deling, the effect of changing the parameters of the object on the accuracy of estimating its coordinates by the observer was investigated. The use of an optimal controller into the automatic control system for the movement of the tunneling shield made it possible to reduce the error in determining the coordinates of the tunne ling shield's frontal point and feedback coefficients when changing the control action, and thereby increase the accuracy of driving the heading board in accordance with the design direction.

Keywords: microtunneling; tunneling shield; observer; motion process; synthesis; optimal regulator; numerical simulation.

Введение

Микротоннелепроходческий комплекс (МТПК) является основным звеном в технологической цепи строительства коммуникационных тоннелей различного назначения. Строительство таких тоннелей ведется в сложных инженерно-геологических условиях с большим диапазоном изменения физико-механических свойств грунтов. Главной задачей при строительстве тоннелей с использованием проходческих комплексов является обеспечение точности выполнения работ в соответствии с проектом.

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

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

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

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

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

Одним из возможных путей оценки недоступных непосредственному измерению координат вектора состояния является применение наблюдающего устройства [2 - 5]. Наблюдающее устройство дает возможность на основании управляющего сигнала и сигнала, пропорционального выходной координате, оценить все остальные координаты вектора состояния.

Синтез наблюдающего устройства

В современных тоннелепроходческих комплексах для передвижения щита применяются гидравлические схемы с объемной гидропередачей, включающей гидравлический насос и гидроцилиндры. Линеаризованные уравнения, описывающие движение щита с гидроприводом на плоскости, имеют следующий вид [1]:

х0 = а01х1 + а02х2>

хх=х2; (1)

х2 = + а22х2 + М^?

в2 в Я

где % = V;а02 = ыаъа21 = -—;а22 =-—; ц=—,

V - линейная скорость движения щита; al - рас-

ISSN 0321-2653 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION.

TECHNICAL SCIENCE. 2019. No 4

" 0 a0i a02 " 0"

где A = 0 0 1 , B = 0

0 a21 a22 _ _y _

стояние между ножевой точкой щита и точкой разворота; х - координаты точек оси щита; Б\ -коэффициент вязкого трения; В2 - коэффициент, учитывающий параметры щита; Л - момент инерции; Я - сила, развиваемая гидроцилиндрами передвижки щита; и - управляющее воздействие; кй - коэффициент пропорциональности.

В матричной форме уравнения принимают следующий вид:

, C = [1 0 0].

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

п

и = X Ьх , где кг - элементы матрицы регулято-

I=1

ра, необходимо иметь информацию обо всех координатах вектора состояния щита. При проектировании системы управления необходимо предварительно оценить его управляемость и наблюдаемость. Для этой цели, на основании системы уравнений (1) могут быть составлены матрицы управляемости (у и наблюдаемости (н и по их структуре выполнена оценка:

Ц(^01 + а02 • а22 ) бу = 0 Ц Ц • а22

Ц( а21 + а22) ] (2) 0

Qh =

0 У • a02

0 y

y У ••

1 0

0 a01

0 a02 (

^02 • l

( a01

^02 ' a22

a22 )

Для оценки недоступных непосредственному измерению координат проходческого щита (угла поворота а и угловой скорости ю) синтезировано наблюдающее устройство [6]. Система уравнений «объект - наблюдающее устройство» имеет следующий вид:

х0 = а0\х\ +а02х2'

х\ ~ХЪ

х2=а2ххх+а22х2 + \1и; (3)

= Х4 + /0^01 (л^ - ) + /о$02 (х2 — Х4);

+ /1^02 (Х2 - Х4 ) + Ци, где хз, Х4 - координаты вектора состояния наблюдающего устройства; /о, /1 - элементы матрицы Ь связи объекта с наблюдающим устройством.

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

.Х1 = Х3; Х2 = Х4 .

Выбором матрицы I = [/0 /1 ]

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

Численное моделирование наблюдающего устройства выполнено в программном комплексе ОрепМойе/1оа. При моделировании исследовано влияние изменения параметров объекта на точность оценки его координат наблюдающим устройством. При этом система уравнений (3) была решена при различных значениях коэффициентов ау. В результате моделирования установлено, что при изменении скорости движения щита V и расстояния от мгновенного центра поворота до ножевой точки щита наблюдающее устройство абсолютно точно оценивает координаты х\(?) и Х2(0.

Рис. 1. Структурная схема объекта и редуцированного наблюдающего устройства второго порядка / Fig. 1. Structural diagram of an object and a reduced second order observer

ISSN 0321-2653 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION.

TECHNICAL SCIENCE. 2019. No 4

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

равной 11,2 % (%), при этом координата X2 (%)

оценивается точно (рис. 2).

% -10"4, рад. ^-10~6,рад.

Wi

Wo

0 20 40 60 80 100 120 140 160 180 200 220 t, с

Рис. 2. Зависимость ошибок оценивания координат и *2 от изменения коэффициента передачи объекта / Fig. 2. The dependence of the estimation errors of the coordinates and *2 on the change in the transmission coefficient of the object

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

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

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

щита микротоннелепроходческого комплекса

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

Условием полной управляемости объекта по критерию Калмана является равенство ранга его матрицы Q1 порядку n объекта. Матрица управления (2) имеет ранг, равный трем. Это указывает на возможность построения регулятора, обеспечивающего любое желаемое располо-

жение корней характеристического полинома системы управления. Для синтеза оптимальной системы автоматического управления воспользуемся принципом Л.С. Понтрягина [7, 8]. Предположим, что оптимальное управление является кусочно-непрерывной функцией |и| < итах .

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

да

I (и) = / (ад2 + 42 п22 + ъ п2 + и2 )№ , (4) 0

где qi (/ = 1, 2, 3) - весовые коэффициенты; П (/ = 1, 2, 3) - координаты возмущенного движения объекта, причем Ц1 = х0 — ; Ц2 = х1 — ; П3 = х2 — х2 х* (/ = 1, 2, 3) - координаты желаемого движения объекта.

Требуется перевести систему из положения Т|1 (0) = Т\1 при / = 0 в положение г|1=0 при

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

Уравнения (1) возмущенного движения объекта запишем в виде

лг =л3;

где al2 = aol, alз = ao2, aз2 = a2l, aзз = a22.

Запишем оптимальный закон управления:

и = 4^1 + ^ + ЬзПз ) ,

где Ц (] = 1, 2, 3) - коэффициенты обратных связей регулятора.

Для определения коэффициентов Ц а^, 2, 3) оптимального закона управления воспользуемся двумя формами записи характеристического уравнения оптимальной системы. Допустим, что корни Хь Х2, Х3 характеристического уравнения известны. Тогда характеристическое уравнение оптимальной системы можно записать в следующем виде:

(p + \ )(p + Х2 )(p + ) = 0.

(5)

После преобразований уравнение (5) может быть представлено в виде

РЪ + У1Р2 + У2Р + Тз = 0, где У1 =Х +Х2 +Х3; у2 +Х2Х3 +ХД2; Уз =к{к2Х3 .

ISSN0321-2653 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION.

TECHNICAL SCIENCE. 2019. No 4

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

р -а12 -а,

det( pI -D1):

0

"43 -1

У?! (y^2 - а32) [Р + (yh, - а33 )]

. (6)

Раскрыв определитель (6), получим характеристическое уравнение

Р3 + P2(yh, -а33) +

(7)

+р(цк2 - а32 + цй1а13) + цк1а11 = 0.

Сравнив коэффициенты при одинаковых степенях р в уравнениях (6) и (7), получим систему алгебраических уравнений:

Ц^з - азз = У^

yh2 - а32 + у?1а13 = у2;

(8)

Ц^1а12 = Уз •

После преобразования из системы (8) получим выражения для определения М, h2, Нз:

У^ _а12 (У2 + а32 ) - Ъа13. _У1 + а33 Ца12 Ца12 Ц

Структурная схема объекта с регулятором, которая справедлива для движения щита как в вертикальной, так и в горизонтальной плоскостях, показана на рис. 3.

h =

Рис. 3. Структурная схема проходческого щита с оптимальным регулятором / Fig. 3. Structural diagram of a tunneling shield's model with an optimal regulator

Численное моделирование процесса движения проходческого щита МТПК

Численное моделирование процесса движения щита [9, 10] проведено при следующих исходных данных: масса щита - 10000 кг, скорость движения - 0,041 м/с, усилие гидроцилиндров передвижки 2000 кН, диаметр щита - 1,42 м, длина щита - 4 м, момент инерции щита относительно оси - 14590 кг-м2.

Численная модель, выполненная в программном комплексе ОрепМойеНеа, представлена на рис. 4.

U

ao2

252

a aoi

t=10 s

Рис. 4. Численная модель проходческого щита с оптимальным регулятором / Fig. 4. Numerical model of a tunneling

shield's with an optimal regulator

Z

w

ISSN 0321-2653 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION.

TECHNICAL SCIENCE. 2019. No 4

В результате моделирования получены графические зависимости изменения угловой скорости а у, угла поворота а и координаты ножевой точки щита Zl при смещении радиус-вектора т2 точки приложения усилия гидроцилиндров передвижки щита. (рис. 5, 6).

а, рад юу, 1/с

0,4

0,2

0,1 0

a

N

Юу

0 10 20 30 40 50 f. с

Рис. 5. График изменения угловой скорости и угла поворота щита при смещении радиус-вектора точки приложения усилия гидроцилиндров передвижки щита / Fig. 5. Dependence of the tunneling shield's angular velocity and angle of rotation on the displacement of the radius vector of the application's point of the hydraulic cylinders force

rz, Zl, мм

0 10 20 30 40 50 t. с

Рис. 6. График изменения координаты ножевой точки щита и смещения радиус-вектора точки приложения

усилия гидроцилиндров передвижки щита / Fig. 6. Changing the coordinate of the tunneling shield's frontal point and the displacement of the radius vector of the application's point of the hydraulic cylinders force

В результате моделирования установлено, что отклонение координаты ножевой точки щита Zi при изменении управляющего воздействия rz составляет 4 мм.

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

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

Вывод

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

Литература

1. Математическая модель пространственного движения проход-

ческого щита микротоннелепроходческого комплекса / А.В. Батюков, А.А. Гуммель, Н.А. Глебов, М.А. Земляной // Изв. вузов. Сев.-Кавк. регион. Техн. науки. 2018. N° 4. С. 72 - 78.

2. Бестраншейная прокладка коммуникаций с применением мик-

рощитов // [Электронный ресурс]. Режим доступа: http:www. sibdom.ru/article.php?id= 122 (дата обращения 01.10.19)

3. Бренер В.А. Щитовые проходческие комплексы: учеб. пособие. М.: 2009. 447 с.

4. Рыбаков А.П. Основы бестраншейных технологий (теория и

практика): Технический учебник-справочник. М.: 2005. 304 с.

5. Глебов Н.А. Управление движением мехатронных комплексов в подземном строительстве: материалы 5-й науч.-техн. конф. «Мехатроника, автоматизация, управление», г. Санкт-Петербург, 14 - 16 октября 2008 г. Спб, 2008. С. 170 - 173

6. Глебов Н.А. Управление движением мехатронных комплексов в неформализованной внешней среде // Изв. вузов. Сев.-Кавк. регион. Техн. науки. 2008. Спец. выпуск; Проблемы мехатроники. С. 58 - 63.

7. Математическая теория оптимальных процессов / Л.С. Понт-рягин, В.Г. Болтянский, Р.В. Гамкрелидзе, Е.Ф. Мищенко. М.: Наука, 1969. 384 с.

8. Воронцов Г.В. Введение в математическую теорию опти-

мального оценивания и управления состоянием технических систем: учеб. пособие. Новочеркасск: ЮРГТУ(НПИ), 2006. 308 с.

9. Рыжиков Ю.И. Решения научно-технических задач на пер-

сональном компьютере. Спб.: КОРОНА принт, 2000. 272 с

10. Fritzson P. Principles of Object-Oriented Modeling and Simulation with Modelica 3.3: A Cyber-Physical Approach 2nd Edition. Wiley_IEEE Press, 2014. 1265 p.

ISSN 0321-2653 ИЗВЕСТИЯ ВУЗОВ. СЕВЕРО-КАВКАЗСКИЙ РЕГИОН._ТЕХНИЧЕСКИЕ НАУКИ. 2019. № 4

ISSN 0321-2653 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION. TECHNICAL SCIENCE. 2019. No 4

References

1. Batyukov A.V., Gummel' A.A., Glebov N.A., Zemlyanoy M.A. Matematicheskaya model' prostranstvennogo dvizheniya prohod-cheskogo shhita mikrotonneleprohodcheskogo kompleksa [The mathematical model of the microtunelling complex's shield's movement]. Izv. vuzov. Sev. - Kavk. region. Tehn. nauki, 2018, no. 4, pp. 72 - 780 (In Russ.)

2. Bestransheynaya prokladka kommunikaciy s primeneniem mikroshhitov [Trenchless laying of communications using microtunnel-ing shields]. Aviable at: http://www.sibdom.ru/article.php?id=122 (accessed 01.10.19)

3. Brener V.A. Shchitovyeprohodcheskie kompleksy: Uchebnoeposobie [Shield tunneling complexes: Tutorial]. Moscow, 2009, 447 p.

4. Rybakov A.P. Osnovy bestransheinykh tekhnologii (teoriya ipraktika): Tekhnicheskii uchebnik-spravochnik [Basics of trenchless technologies (theory and practice): Technical reference manual]. Moscow: Publ. PressByuro, 2005, 304 p.

5. Glebov N.A. [Motion control of mechatronic complexes in underground construction]. Materialy 5-i nauchno-technicheskoi konferentsii «Mekhatronika, automatizatsiya, upravlenie» [Materials of the 5th scientific and technical conference "Mechatronics, automation, control"]. Sankt- Petersburg, 2008, pp. 170 - 173. (In Russ.)

6. Glebov N.A. Upravleniye dvizheniyem mekhatronnykh kompleksov v neformalizovannoy vneshney srede [Motion control of mechatronic complexes in an informal external environment]. Izv. vuzov. Sev. - Kavk. Region. Tech. nauki, 2008, Special issue. The problems of mechatronics, pp. 58 - 63. (In Russ.)

7. Pontryagin L.S., Boltyansky V.G., Gamkrelidze R.V., Mishchenko E.F. Matematicheskaya teoriya optimal'nykh protsessov [The mathematical theory of optimal processes]. Moscow: Nauka, 1969, 384 p.

8. Vorontsov G.V. Vvedeniye v matematicheskuyu teoriyu optimal'nogo otsenivaniya i upravleniya sostoyaniyem tekhnicheskikh sistem: Uchebnoye posobiye [An Introduction to the Mathematical Theory of Optimal Estimation and Management of the Condition of Technical Systems: Textbook]. Novocherkassk: SRSTU (NPI), 2006, 308 p.

9. Ryzhikov Yu.I. Resheniya nauchno-tekhnicheskikh zadach na personal'nom komp'yutere [Solutions of scientific and technical problems on a personal computer]. Sankt-Petersburg: KORONA print, 2000, 272 p.

10. Fritzson P. Principles of Object-Oriented Modeling and Simulation with Modelica 3.3: A Cyber-Physical Approach 2nd Edition. - Wiley_IEEE Press, 2014, 1265 p.

Поступила в редакцию /Received 13 ноября 2019 г. /November 13, 2019