Моделирование динамики систем с переменной массой и управление производственными системами Текст научной статьи по специальности «Математика»

УДК519.711.3, 531, 531.391.5

Р. Г. Мухарлямов, Н. В. Абрамов, А. А. Ахметов

МОДЕЛИРОВАНИЕ ДИНАМИКИ СИСТЕМ С ПЕРЕМЕННОЙ МАССОЙ

И УПРАВЛЕНИЕ ПРОИЗВОДСТВЕННЫМИ СИСТЕМАМИ

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

предприятие.

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

Keywords: dynamics, control, equation, constraints, stabilization, numerical solution, system.

Methods of mathematical modeling of dynamic objects allow you to apply the general positions of classical mechanics to study the dynamics of physical systems and control problems of industrial enterprises. For describer the dynamic processes in economic and production systems equations of analytical system dynamics with variable mass used. The conditions of stability of dynamics equations solutions of the manufacturing systems with respect to the constraint equations that reflect control objectives, and the algorithm of numerical solution of control are defined.. The problem ofplanning and control of industrial enterprise for the production of tires is considered.

Введение

Методы динамики и математического моделирования [1] существенно расширили область применения общих положений классической механики для исследования систем различной природы и решения задач управления производственными предприятиями [2],[3]. Современные динамические аналогии [4],[5],[6],[7] позволяют использовать уравнения классической механики для описания динамических процессов управляемых систем с соблюдением уравнений связей. Так для описания динамических процессов в экономических объектах и производственных системах могут быть использованы уравнения аналитической динамики систем с переменной массой [8]. Необходимым условием обеспечения связей является устойчивость решений уравнений динамики по отношению к уравнениям связей. Для стабилизации связей при численном решении уравнений динамики [9] на управляющие воздействия должны быть наложены дополнительные ограничения, определяемые выбором метода численного решения [10].

Моделирование динамики управляемой системы с переменной массой

Динамика управляемой системы с переменой массой описывается дифференциально-алгебраическими уравнениями динамики несвободной механической системы [8]

ч =—,

сИ

d д L д L

----=0+Ф,

dt 5q 5q

д^д) = 0, д (ц, Ч, ^ = 0, (2)

чОо) = ч0, Ч(^) = Ч° . (3)

Здесь д = ..,%) , Ч = (1 ..,С1П)

— векторы обобщенных координат и скоростей; I ) = Т-Р - функция

Лагранжа, 2Т = Чт А^, ^4 +2а т(д, 1) <Ч + а0 (^,о) -

удвоенная кинетическая энергия, П = п(|)

потенциальная энергии системы, О, Ф — векторы непотенциальных и реактивных сил,

д(сд) = (д1,...,дт), д(ч^д)(дт+1,^,дг) —

левые части векторов голономных и неголономных связей, наложенных на систему. Соотношения (3) задают начальные условия.

Допуская возможные отклонения от уравнений голономных связей, их производных по времени и уравнений неголономных связей при численном решении системы (1)-(3), заменим уравнения (2) уравнениями программных связей

ш(д, д, ^ = г,

(4)

w =

" дМ g (q, q, t)

g(q,q, t)

f

g = Gq + gt,

M =

sg,

л

_M

V

at

j

= 1,...,т , ! = 1 ,...,п .

Под невозмущенным движением системы будем понимать движение вдоль многообразия,

определяемого уравнениями связей (3). Любые другие движения будут возмущенными.

Вектор ъ = (ъ.,..., ът+г) оценивает

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

dz

— = Pz, dt

Здесь P =

vj^ I P10 P11 P12.

(5)

- квадратная матрица, I

- единичная матрица (т х т); 01, О2 - нуль-матрицы ^ш), (т х(г - т)) ; Р10, Р., 1 -вещественные матрицы размерности (г х т), Р.2 —

матрица (г х (г - т)).

Предполагая, что программные связи (4) являются идеальными, представим уравнения (1)-(3) системой дифференциальных уравнений [10]

— = а,

С11

ас1 -1/ т ~Т \

— = А (F + G Л + G Л),

dt v '

q(t0) = q0, q(t0) = q°.

(6)

где G = (gpi) =

f Л v5 q У

r = m + 1,...,r; А, Л -

векторы множителей Лагранжа;

д

Ф = 0 + Ф-Ад-а + —.

Дифференцирование уравнений (2) по времени с учетом (5) позволяет определить вектор множителей Лагранжа:

-

Л II _D1 d2"

Л _D3 D4_

(P1Z-F)

(7)

P1 [P10 P11 P12 ]

D1 = GA-1GT, D2 = GA-1GT,

~ -1 T ~ -1 ~ T

D3 = G AG, D4 = G AG,

Ф

GA'^0 + (Gq) q + 2<3tq + g

GA Ф + Gq + gt

f

n 5 9„ .

=—q

Vl=у

С 2 Л

Г \

5gp

V^i У

f 2 Л

8 gm

V Л У

Л( =

f л \

здр

V 5t У

Правые части системы уравнений (6) при

определенных выражениях Л, Л представляют

функции переменных q,q. Окончательный вид

правых частей уравнений (6) определяется выбором

матриц Рю>рц>р12, которые можно использовать

для обеспечения устойчивости решения системы относительно уравнений связей (2).

Построим вычислительную схему решения задачи (6). Будем аппроксимировать непрерывно

дифференцируемые векторы q = q(t), ] = ](1),

кусочно-непрерывными (ступенчатыми) функциями

С (1) = С (1к ) ,

чОНК)■

к = 0.....N -1, ^ = Т.

При заданных начальных условиях С ( ^ ) = ]0, с] ( ^ ) = с]0 итерационный процесс

решения системы уравнений (6) сводится к построению разностных уравнений метода Эйлера " (к+1) (к) .(к)

.(к+1) .(к) ..(к) I а4 - qw + <ач т,

(8)

1к+= kт=

т = const.

Устойчивость решения

системы (6) относительно уравнений связей (2) определяется устойчивостью тривиального решения системы (5).

Используя метод функций Ляпунова У = У (ъ),

сформулируем условие устойчивости численного решения уравнений динамики (6).

(к)

Обозначим V — значение функции Ляпунова V = V (ъ) в момент времени 1 = 1к . определим условия соблюдения неравенства

у(к+1) < у(Ю . (9)

Положим 2V = ъТКъ, где К - постоянная

(к+1) (к) (к) (к) матрица, и . Величина

определяется из соотношений (6),(6),(8):

AV(k)=(zTKPz)(k)T + ±(zTPTKPz)'

(Ю 2

т + .

Считая величину Т достаточно малой,

(к)

условие выполнения неравенства представим в упрощенном виде:

( T \(k)

i z KPz) <0, k = 0,...,N. (10)

Выполнение условия (10) может может быть обеспечено соответствующим выбором коэффициентов матрицы P^ при фиксированном шаге интегрирования т. Таким образом, имеем задачу определения такой матрицы P.|, чтобы ограничение (10) выполнялось во всех узловых точках tk, k = 0,...,N.

Эту задачу можно решить либо начальным выбором соответствующих значений P| (tk ) , что

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

Представим левую часть неравенства (10) в виде функционала

|[р] = ma (zTKPzfk). L 1J ke[0,K]V '

Введем меру близости р[Р]] до области

Q : I J < 0 . Мера р характеризует удаление

механической системы от множества допустимых решений

рЬН

[о, 1<0.

Тогда условие (10) примет вид

р[р^] = 0. (11)

Таким образом, найденное допустимое решение и соответствующие параметры системы

обеспечивают минимум функционала р [р^ J.

Алгоритм численного решения. Порядок численного решения системы (6),(6), при ограничении (11) можно представить блок-схемой, представленной на рис. 1.

Рис. 1 - Блок-схема алгоритма численного решения уравнений динамики

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

Вектор ДР/1

точки

P[Р] к точке P[Р+1]

1 Л IVUW I 1

направлением антиградиента

ДрИ= -£Иу|[р],

направления перехода от можно представить

где

[р]

Г ЛИ

51

д У

- градиент функционала

р[р] в точке P1 = P1[p], P1[p] к точке p/P+1] .

.И

- шаг перехода от

Моделирование динамики простого производственного объекта

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

которого выразим величиной х (1) . Под потоком

продукции будем понимать объем произведенной

продукции в единицу времени ( \

х ( ) =

произведенной

dt

где X (1) - объем

продукции к моменту времени 1.

Величину X (1) будем измерять в

количественном выражении - шт., тн. и т.д. Тогда

величина х (1) будет иметь размерность шт./год,

шт./мес. и т.д.

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

X = X - а (1) X,

где а - коэффициент затрат на рецикл. Пусть простой производственный объект обладает основными производственными фондами

(ОПФ) в количестве 8(ф) () , ф = 1.....ф и

снабжается извне оборотными фондами (ОФ)

0 ,п=1.....п

Компоненты вектора

/ (1) (2) (Ф)\

ОПФ - средства труда,

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

первоначальной стоимости без учета амортизации.

( 0) (2) (П)\ Компонентами вектора V = I V .V ,...,м I ОФ

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

Компоненты вектора V также могут измеряться в количественном или стоимостном выражении в единицу времени. Векторы ОПФ и ОФ соответствуют невозмущенному развитию производственного объекта, если они являются векторами с пропорциональными коэффициентами [11]

,(2) _(Ф)

3 у(),

т

к(2)

>)

т

„(п)

(5)

кп) к1)

где т Ф (1), к (^ — мгновенная фондоемкость и материалоемкость продукции объекта по ОПФ

вида ф и ОФ вида тт соответственно, к( ) -трудоемкость продукции, Ь - численность основных рабочих; у — производственная мощность объекта

у () =—, dt

где У (1) - максимально возможный выпуск

продукции к моменту времени 1 при нестесненных ОФ объекта. Производственная мощность измеряется в тех же величинах, что и поток выпуска продукции. Справедливо соотношение [6]

0^х(1:)^у(1:), означающее, что поток выпуска

продукции не превышает производственной мощности объекта.

(ср) (П) (Ь)

Коэффициенты , , к означают

эффективность использования ОПФ, ОФ и людских

ресурсов. Уменьшение значений коэффициентов

т(ср^, к(тт), к(Ь) влечет более интенсивное использование ОПФ, ОФ и труда, когда увеличивается количество продукции, получаемое с каждой единицы физического капитала и оборотных средств. Увеличение объема выпуска продукции за счет улучшения использования ОПФ и ОФ равносильно расширению производства (приросту мощностей) без дополнительных капитальных вложений на строительство новых

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

Рассмотрим производственный объект по производству автопокрышек мощностью у' = 400000 шт. / год (у' - поток товарной продукции). В качестве основного технологического оборудования выступают сборочные станки в

количестве

0)

,

окрасочные станки - э

(2)

форматоры-вулканизаторы - э

,(4)

(3)

обрезки выпрессовок - э

станки для и станки для

„( 5)

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

Таблица 1 - Характеристики ОПФ

Наименовани Коэффициен Производи Эффективн

е ОПФ т возвратных а (Ф) потерь, тельность единицы ОПФ, Г шт/час ый годовой фонд времени, Т, час

1 2 3 4

Сборочный 0,0015 2,4 7781

станок

Окрасочный 0,0011 393,8 7872

станок

Форматор- 0,0011 0,667 7701

вулканизатор

Станок 0,0001 40,8 7890

обрезки

выпрессовок

Балансировоч 0,00005 43 7656

ныи станок

С учетом суммарного коэффициента возвратных потерь

ф (ф) а = 2 а^ = 0,00385

Ф=1

рассчитаем фондоемкость продукции по видам ОПФ и состояние ОПФ по формулам

(ф) 1 (Ф)

> = —, >

Т1

Результаты сведем в табл. 2.

т

0 + а)

Таблица 2 - Расчет основного оборудования

Наименование Коэффициент Состояние

ОПФ фондоемкости, гг/4^, ед млн.шт/ год ОПФ, , ед

Сборочный 53,6 21,5

станок

Окрасочный 0,323 0,13

станок

Форматор- 194,7 78,2

вулканизатор

Станок 3,11 1,25

обрезки

выпрессовок

Балансировочн 3,04 1,22

ый станок

Зная первоначальную стоимость каждого станка можно выразить состояние ОПФ в денежном

руб

выражении, а фондоемкости в---.

шт/ год

Пусть при производстве автопокрышек

0) (2)

необходимы следующие виды ОФ , -

резиновые смеси двух видов, -

(4) й (5)

электроэнергия, - сжатый воздух, -

(6)

оборотная вода, V - перегретая вода, \ -

численность рабочих.

Пусть даны нормы расхода ОФ

производства и норма обслуживания. Потребные

ОФ определим по формулам

(тт) (тт) (Ц

, .

Результаты сведем в табл. 3.

Таблица 3 - Расчет ресурсов (потока материалов и численности)

Наименование ОФ Нормы расхода, к« Поток ОФ V1 ' и численность Ь,

1 2 3

Резиновая смесь I, 2 / 26,6 м /шт 1219 м2/час

Резиновая смесь II, 48,5 кг/шт 2223 кг/час

Электроэнерги я, 1,15 кВт/шт 52,7 кВт/час

Сжатый воздух, 5,57 м3/шт 3 / 255 м /час

Оборотная вода, 0,823 м3/шт 37,7 м3/час

Перегретая 0,77 м3/шт 35,3 м3/час

вода,

Трудозатраты, чел. 0,428 тыс.шт/ год 171 чел

Если компоненты вектора ОПФ измеряются в условных денежных единицах, а вектор потока ОФ - в условных единицах в единицу времени, то

величины векторов Б (1), У^) производственного

объекта можно измерить в виде суммы своих компонентов [6]

3 (Ф) И (тт)

э = 2 ', V = 2 V1 '. ф=1 тт=1

Состояние ОПФ в стоимостном или количественном выражении изменяется во времени в связи с выбытием оборудования без учета износа, реализации другим производственным объектам, а также капитальным вложениям. Предполагая, что скорость выбытия пропорциональна состоянию ОПФ [6]

— = с11

где

Б1 —

состояние ОПФ,

эксплуатации, уравнение развития производственного объекта имеет вид

ds ds1 / \ - + —= и1 ()

выбывших из ОПФ

dt

или

dt

ds

— = -ьб., + dt

и ()

(12)

выбытия ОПФ (ликвидации т.д.), и1 (1) -

где

- фондоемкость продукции

где - коэффициент

производственного объекта

оборудования из-за износа и

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

и1 (1) > 0 либо выбытие ОПФ в результате

передачи (реализации) оборудования другим

производственным объектам и (1) <0 .

Состояние в ОПФ определяет производственную мощность объекта [12].

б (1) = т (1) у (1), (13)

ф (<р) 2 т1^

Ф=1

производственного объекта по всем видам ОПФ.

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

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

т (1) = т0+т1 (1 )-т2 (),

где то - начальная фондоемкость; т1 - величина

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

отражаемые в состоянии Б; т2 - величина

фондоемкости, которая убыла к моменту 1, вследствие внедрения научно-технического прогресса, совершенствования технологии производства и т.д.

С учетом (12), (13) уравнение динамики ОПФ примет вид уравнения развития мощности производственного объекта

с1т0

с1у (Зт,,

т— = -(Вту + и,,---у + —— у.

с11 с11 с11

Полученное уравнение совпадает

(14)

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

Р = -

dm.

dm-

У +

dt dt

есть аналог реактивной переменных масс.

силы

механике

с

в

Аналогия развития мощности

производственного объекта с динамикой материальной точки переменной массы позволяет применить методы аналитической динамики в управлении производством. Тем не менее, эта аналогия чисто математическая [11] и в случае с динамикой производственных объектов слагаемые, составляющие величину Я, требуют экономического толкования.

dm,

Слагаемое

'1

dt

у измеряется в единицах

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

Если скорость накопления фондоемкости

dm„

задана, то выражение

1

dt

у можно отнести к

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

с1т.

'1

сИ

у = уу():

где - коэффициент морального и физического износа.

dm2

Величина -у характеризует экономию

dt

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

dm.■

'2

dt

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

в совершенствование эксплуатации оборудования dm.■

2

dt

= и2 (t).

(15)

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

1. инвестиции, обеспечивающие восполнение и наращивание мощностей при неизменной

технологии производства - и ( ) ;

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

и2(t).

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

(1-а) У = Р( t)

(16)

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

Вместо цели управления (16) будем использовать уравнение программной связи

(1- а)у-р(^) = а(1), (17)

где правая часть уравнения (17) будет определяться из системы уравнений возмущений связей

— = а, сИ

dа

А( У,а,аД),

где А(у,У,0,0Д) = 0 .

В качестве А (у, У, а, а, 1) примем

линейную функцию

С - с (у,Уд)а + с2(у,У,1)а.

При построении математической модели развития производственного объекта потребуем,

чтобы управляющие воздействия ы1, ^

обеспечивали устойчивое развитие объекта относительно цели управления.

Параметр С1, С2 должны быть выбраны

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

и - ш(у,У,а,ад) + —, 2

и- - -ш(у,У,а,ад) + —, 2

где со = и)(у,У, 1) - произвольная функция,

-(у + Уа + р + с,,а + с2 ((1 - а)у - '/а -р))

т

(а)

+(3ту + уУ

Работа выполнена при финансовой поддержке РФФИ, номер проекта 13-08-00535.

Литература

1. Мухарлямов Р.Г. Моделирование динамических процессов различной природы // Проблемы аналитической механики и теории устойчивости. Сборник трудов, посвященный памяти академика В.В. Румянцева. Изд-во «Наука».. 2009. С. 310-324.

2. Сиразетдинов Т.К. Динамическое моделирование экономических объектов. - Казань. - «Фэн». - 1996. -223 с.

3. Сиразетдинов Т.К., Родионов В.В., Сиразетдинов Р.Т. Динамические модели экономического региона. -Казань: Изд-во «Фэн» Академии Наук РТ, 2005. - 320 с.

4. Ольсон Г. Динамические аналогии. М.: Гос. изд-во ин. лит-ры, 1947 - 224 с.

5. Самарский, А. А. Математическое моделирование. Идеи. Методы. Примеры. Изд. 2-е, исправл. — М.: Физматлит, 2002. - 320 с.

6. Layton R. Algebraic-Differential Equations of Mechanical Systems / Springer. 2001. 159 p.

7. Мухарлямов Р.Г. Стабилизация движений механических систем на заданных многообразиях фазового пространства // ПММ. 2006. Том 70. № 2. С 236-249.

8. Н.В. Абрамов, Р.Г. Мухарлямов. Ж.К. Киргизбаев. Управление динамикой систем с программными связями: Монография. - Нижневартовск. Изд-во Нижневарт. гос. ун-та, 2013. 202 с.

9. Baumgarte J, Stabilization of constraints and integrals of motion in dynamical systems // Comp. Math. Appl. Mech. Eng. 1 (1972), 1-16.

10. R.G. Mukharlyamov, Assaye Walelgn Beshaw. Solving Differential Equation of Motion for Constrained Mechanical Systems // Вестник Российского университета дружбы народов. Серия «Математика. Информатика. Физика». № 3, 2013. С. 81-92.

11. Р.Г. Мухарлямов, О.В. Матухина. Моделирование процессов управления, устойчивость и стабилизация // Вестник Казан. технол. ун-та, 2012, т. 15, в. 12. - С. 220-224.

12. Л. Н. Дюкалов, Ю. Н. Иванов, В. В. Токарев. Теория управления и экономические системы: I // Автоматика и телемеханика. — 1974. — № 5. — С. 117-132.

13. Матухина О. В. О моделировании динамики в информационных системах символьных вычислений // Вестник Казан. технол. ун-та, 2013, т. 16, в. 2. - С. 191-193.

© Н. В. Абрамов - кафедра естественнонаучных дисциплин Филиала Тюменского Государственного нефтегазового университета [email protected]; А. А. Ахметов - старший преподаватель кафедры математики НХТИ КНИТУ; Р. Г. Мухарлямов - профессор кафедры математики НХТИ КНИТУ, [email protected].

© N. V. Abramov - chair of natural-science disciplines of Branch of Tyumen State Oil and Gas University of [email protected]; A. A. Akhmetov - the senior teacher of department of mathematics of КСНТ1 KNRTU; R. G. Mukharlyamov - professor of department of mathematics of ЫСНТ1 KNRTU, [email protected].