CC BY
19
--------------------□ □------------------------
Описана ідея спільного застосування методу топологічної згортки та оптиміза-ційного розрахункового алгоритму для розрахунку гідравлічних мереж, які містять висячі вузли. Наведені теоретичні основи методів, їх основні переваги, а також практичний приклад розрахунку за допомогою запропонованої методики
Ключові слова: гідравлічна мережа, висячий вузол, оптимізація, розрахунок
□--------------------------------------□
Описана идея совместного использования метода топологической свертки и оптимизационного расчетного алгоритма для расчета гидравлических сетей, содержащих висящие узлы. Приведены теоретические обоснования методов, их основные достоинства, а также практический пример расчета с помощью предложенной методики Ключевые слова: гидравлическая сеть, висящий узел, оптимизация, расчет
□--------------------------------------□
This article represents the idea of joint use of such methods as the method of topological furl and the optimization calculation algorithm for hydraulic networks’ with dangling nodes parameters calculation. The theoretical basis of mentioned methods, as well as the practical example of calculation by means of the offered methodology, are given
Keywords: hydraulic network, dangling node, optimization, calculation --------------------□ □------------------------
УДК 519.853
ОПТИМИЗАЦИЯ РАСЧЕТА ГИДРАВЛИЧЕСКИХ СЕТЕЙ С ВИСЯЩИМИ
УЗЛАМИ
С.А. Исаенко
Магистрант* Контактный тел.: 067-456-00-07 Е-mail: svetlyachok7@gmail.com
В.Н. Медведева
Кандидат технических наук, доцент* Контактный тел.: (044) 400-85-45 Е-mail: romyald@voliacable.com
Ю. Д. Ще рбашин
Кандидат технических наук, доцент* Контактный тел.: (044) 454-95-23 Е-mail: romyald@voliacable.com *Кафедра автоматизации проектирования энергетических
процессов и систем Национальный технический университет Украины «Киевский политехнический институт»
Введение
В процессе проектирования и реконструкции гидравлических сетей или в случае возникновения аварийных ситуаций часто возникает острая необходимость оперативных численных расчетов гидравлических сетей для правильного выбора новых конфигураций, подбора и размещения оборудования таким образом, чтобы при выполнении требуемых технологических ограничений экономические потери были минимальными. Поэтому основной целью современных исследований в этой области должно быть повышение надежности, эффективности и качества расчета за счет математической оптимизации расчета и настройки гидравлических систем компьютерными средствами.
В данной статье рассматривается идея совместного применения усовершенствованного алгоритма решения задач нелинейного программирования с использованием базиса сокращенной размерности [1] и метода топологической свертки [2] для расчета гидравлических сетей, содержащих висящие узлы. Использование метода топологической свертки позволит решить проблему висящих узлов, часто возникающую при расчете больших магистральных сетей, нефтепроводов, систем пароснаб-
жения, а благодаря усовершенствованному алгоритму решения задач нелинейного программирования с использованием базиса сокращенной размерности, станет возможным произвести расчет сети с использованием математической оптимизации, что, безусловно, повлияет на скорость и эффективность расчета.
Постановка задачи
Говоря о расчете потокораспределения в гидравлических сетях, получим следующую оптимизационную задачу:
Найти
] 0"- (Й ь^Х )2 Й ЙЙ тах, xj, WЙ Еп (1)
]=1
если область допустимых решений О определяется ограничениями
У. =£а^ + 0,1 = 1..(т -1), (2)
]=1
z1 = Ф,(Х) = Щ = 0,1 = 1..(с -1), (3)
]=1
vj = xj - X“in > 0,j = 1..n,
J J J J
vj = —xj + X”“ > 0, j = 1..n,
J j j j
(4)
(5)
где
Хj - элемент вектора неизвестных переменных -расходы на ветвях схемы;
Ф0 - функция цели, которая максимизируется (сформулированная таким образом, что ее максимальное значение равно нулю);
Ц - элементы матрицы контуров системы;
Sj - гидравлическое сопротивление ]-ой ветви схемы;
у. - отклонение в точке X от 1-ой ограничивающей гиперплоскости;
a1j - элемент 1-ой строки матрицы А соединений схемы;
- свободный член - элемент вектора расходов в узлах схемы;
z1 - отклонение в точке X от 1-ой ограничивающей поверхности Ф1(Х) = 0;
Ф1(Х)- 1-тое нелинейное ограничение-равенство;
Хт1П(Хтах) - нижний (верхний) предел допустимого изменения значения неизвестной переменной Хj;
у^) - отклонения в точке X от нижнего (верхнего) предела допустимых значений переменной Хj.
Значения Хт1П(Хтах) принимаем равными очень большому отрицательному (положительному) числу, то есть Хт1п ^-те (Хтах ^+~) . Выбранные таким образом границы гарантируют, что в процессе решения мы никогда не выйдем за эти границы. И, следовательно, нет необходимости в вычислении шага до ограничивающей гиперплоскости тривиальных ограничений.
В качестве начальных приближений для вектора неизвестных переменных нужно выбирать нулевые значения.
Кроме того, в нашей задаче мы рассматриваем гидравлические сети, содержащие висящие узлы, то есть сети вида:
Рис. 1.
Сворачивание висящих узлов
Метод позволяет исключить из расчета висящие узлы графа гидравлической сети.
Для рассмотренной выше топологии гидравлической сети введем обозначения расходов в ветвях и давлений в узлах схемы:
Метод топологической свертки впервые предложен в работе [2] и назван методом линейной свертки. В последующем используется термин топологическая свертка (ТпС).
Р; - давления в узлах; Х; - расходы в ветвях;
Пусть Р - Р; = а; + Ь; ■ Х; (1=0 -^к), где а; (1=0 -^к) - некоторый перепад давления, а Ь; (;=0 ^к) - коэффициент сопротивления ветви.
Зададимся некоторым значением давления Р=Р1 в висящем узле и значением Р1+1. Тогда 1
Хц = Ц (Р1 - Р- а;) (1=1+к);
к 1 к 1 к р + а Х0,=-^(р1-Ч-.,) = -Р1^11|Г+1^ (•>
1
Х12 = — ■ (Р1 +1 - р - а;) (;=1-к);
к ! к ! k р + a
;=—£ 1 .(P1 +1 — р_ — a,) = -(P1 +1).£ 1 + £^ (7)
а для изменения расходов в ветвях при изменении значений давления в висящем узле на единицу получим:
к 1 1
dxo = Х0,2- Х0,1 = ; dxl = Х1,2- Хи=ц (;=1-к);
t-ib,
dx,/dx0 = —y,/£yn (i=1-k)
(8)
Пусть Р0 - значение давления, при котором Х0(Р0)=0. Тогда расход в любой из ветвей системы можно определить как линейную функцию от давления в узле Р:
Х;(Р) = у; (Р0-р -а;) + (Р-Р0)у = к
= Xi,0 — X0(P)
у,/ £ Уп
(i=1 - к)
(9)
Xo(P) = (P — P0). £ yn; P0 = 1 £ у, (P,+ a,) / £ Уп (10)
п=1 V ,=1 ) п=1
Из (10) тогда следует, что для решения задачи расчета расходов достаточно найти расход X0, а поэтому вместо всех ветвей узла достаточно оставить только одну ветвь, по параметрам которой можно определить расход X0.
X
n=1
n=1
Е
Таким образом, мы исключили из расчета висящие узлы схемы. Далее проведем непосредственный расчет гидравлической сети с помощью усовершенствованного алгоритма решения задач нелинейного программирования.
Расчет потокораспределения
Усовершенствованный алгоритм решения задач нелинейного программирования с использованием базиса сокращенной размерности относится к классу методов аппроксимирующего программирования [9],[6] и подробно описан в работе [1]. Как и в линейном программировании, здесь на каждой итерации определяется набор из п эффективных (имеющих нулевые уклонения) ограничивающих плоскостей, образующих базисный конус. Однако, в отличие от симплексных методов (движения по одному из ребер базисного конуса), метод допускает также движение ’’внутрь” многогранника за счет одновременного изменения нескольких уклонений.
Геометрическая интерпретация процесса решения задачи - это движение изображающей точки X® = =со1{Х1(к),...,Хп(к)} в п-мерном евклидовом пространстве Еп:
X(k+1) = X® +ДX(k) = X® + Х(к)о(к) , (11)
где:
о(к) - луч, исходящий из X®, движение вдоль которого не нарушает немедленно удовлетворенные в X® равенства,
Х(к)>0 - длина шага вдоль о(к).
По сравнению с X®, точка X(k+1) должна либо уменьшить (ликвидировать) невязку одного из неудовлетворенных в X® равенств при сохранении ранее удовлетворенных равенств, либо (при нахождении X® в области допустимых решений О) улучшить функцию цели (1).
Итеративный алгоритм (11) требует решения на каждой итерации двух проблем:
1) Выбор подходящего направления - луча о®, вдоль которого уменьшается невязка неудовлетворенного в X® равенства из числа (2),(3) или улучшается функция цели.
2) Определение допустимого шага Х(к), соответствующего одному из условий:
- достигнута гиперповерхность р-го нарушенного в X® равенства Фр^) =0 ,
- нарушена инженерная точность ±£з удовлетворенного в Х(к) равенства Фа^)^,
- достигнут экстремум функции цели вдоль луча о(к).
Нелинейные ограничения Ф;^) в точке X® заменяются их линейной аппроксимацией:
Ф^)^ = а^ЬДК +ф1(X(k)), (12)
где:
ДX=X - X(k)=co1{Дxj, ^=1,..,п} - вектор отклонений от точки X®,
а^^да | х-х(к)={д Ф^/ЭХ) I х=х(к), ]=1,..,п}-градиент Ф;(Х) в точке Х(к) (13)
Градиенты (нормали) плоскостей общего положения (2) и плоскостей тривиальных ограничений-неравенств (4),(5), параллельных координатным плоскостям, постоянны в процессе решения задачи и не требуют дополнительных вычислений типа (13):
Vy1(X)=a1=const•; Уу|' (X)= ej=const•;
Vуj"(X)= (14)
где:
a1 =(а;1,..,а;|,..,а1п)- нормаль ;-ой гиперплоскости (2), в] =(0,0,..1,0..0)- п-мерная вектор-строка с ]-ым единичным компонентом и остальными нулевыми компонентами.
В качестве начального приближения искомых переменных примем их нижнюю границу: X(0)=co1{xj(0), ]=1..п}=то1(Х1тт,..,Хптт), а в качестве базисных переменных - уклонения У]', ] = 1,..,п.
Решение разобьем на 2 этапа:
Этап I - последовательное введение в базис уклонений у;, ;=1,..,т и удержание их в базисе до конца решения задачи;
Этап II - аппроксимация нового эффективного в текущей точке нелинейного ограничения-равенства Фр^®)^; введение в базис Zp вместо zs, для которого Ф^) в текущей точке X(k) стало ненулевым (Iфs(X(k))|>es), или вместо у/(у"), по которому уклонялись.
Контроль уклонений от эффективных ограничений
Для определения о(к) используем в качестве базиса независимых переменных уклонения от эффективных ограничивающих плоскостей [2].
Исходная задача (1)-(5) в качестве базиса имеет координаты Х1,..,Хп искомых переменных, а в качестве зависимых от базиса величин - уклонения у^у^У]" Задача состоит в алгебраическом преобразовании, где зависимыми переменными становятся Х^..,Хп , а независимыми (базисными) - п уклонений из перечисленных выше:
X=G(k)* Y Б(k)+b(k) (15)
Здесь:
X=co1{xj,j=1..n} - вектор искомых переменных,
G(k) - обратная матрица размерности п*п зависимости искомых переменных от уклонений от ограничивающих плоскостей текущего базисного конуса,
YБ(k)=co1{y1, z1, у]', у]"; i,jeБ(k)} - п-мерный вектор уклонений от плоскостей текущего базиса Б(к),
Б(к) - множество индексов базисных плоскостей, dim(Б(k))=n,
b(k) = X« - вектор свободных членов - вершина, образованная пересечением п плоскостей базисного конуса.
Далее в целях упрощения нотации верхний индекс номера итерации опускаем.
В соответствии с (2)-(5) уклонения от базисных При окончании Этапа I обратная матрица имеет
плоскостей зависят от искомых переменных: структуру:
YF
где:
*X +Ьб
(16)
YБ - вектор уклонений от базисных плоскостей; в начальный момент YБ =со1{У)', ]=1,...,п},
Аб - п*п-матрица нормалей плоскостей базисного конуса; в начальный момент Аб=1п - единичная п*п матрица,
X=col{xj, ^=1,.,п} - вектор искомых переменных,
Ьб - вектор уклонения начала координат от плоскостей базисного конуса; в начальный момент Ьб =-Xmm =
=со1{- Х^іп, ]=1,_,п}.
Из (16) нетрудно получить обратную зависимость:
X= G* YБ +Ь, где:
G=(AБ)-1
(17)
; і=1,..,п; j=1,..,n}-
атная п*п-ма-в начальный момент G=In
трица с элементами - единичная п*п матрица,
Ь = - G*bБ
Ь - п-мерный вектор-столбец свободных членов; в начальный момент Ь = Хтіп. Не вошедшие в базис уклонения вычисляются через базисные уклонения подстановкой (17) в (2):
Уі=«і* YБ + V
где:
ui= ai*G Ьі'=Ьі + ai* Ь
(18)
(19)
(20)
На Этапе I последовательно в цикле вводим в базис уь_,уп. Для этого выполняем следующие действия: Шаг 0. Формируем начальную ‘расширенную’ (дополненную вектором Ь) обратную матрицу G/ = ^,Ь)= =(1п,Ь) размерности пх(п+1). Номер цикла t=1.
Шаг 1. По формуле (19) вычисляем вектор-строку ^ и ищем в ней максимальный по модулю элемент иіз. По формуле (20) вычисляем свободный член Ь/ и формируем ‘расширенную’ вектор-строку ^ Ь/).
Шаг 2. Выполняя над расширенной обратной матрицей шаг жорданова преоб-разования [2], вводим в базис у4 , исключая из базиса V/:
V &
&І
Ці = (1,..,n),j = (1,..,П + Ш * 8.
gls = ^,і = 1,..,п,
и(.
где:
(21)
(22)
g1j - элементы ‘расширенной’ пх(п+1) матрицы G', g1j' - элементы скорректированной ‘расширенной’ пх(п+1) матрицы G',
utj - элементы ‘расширенной’ вектор-строки Шаг 3. 1=1+1. Если 1>т, то выход из цикла, иначе
- переход к Шагу 1.
G=(L,H,b),
а вектор искомых переменных X= L*Y + H*V + Ь,
(23)
где:
Y = со1(уі, у2,.., ут) - т-мерный вектор-столбец уклонений от плоскостей-линейных ограничений (2),
V = со1{ Vj/,vp"•, j,pєБ} - (п-т)-мерный вектор уклонений от тривиальных плоскостей, входящих в текущий базис,
Б - множество индексов базисных плоскостей, dim(Б)=n,
Ь - вектор свободных членов - вершина, образованная пересечением п плоскостей базисного конуса,
L,H - подматрицы обратной матрицы G размерности пхт и пх(п-т) соответственно.
Из (23) следует, что в момент окончания Этапа I изображающая точка X(m) - вершина базисного конуса, равна Ь. Используем X(0)=b как начальную точку для Этапа II, где вычисления будем производить в приращениях в соответствии с алгоритмом (11).
Поскольку на Этапе II линейные ограничения-равенства (2) никогда не выходят из базиса, т.е. Y=0 до конца решения задачи, то в дальнейшем можно оперировать только уклонениями V и их приращениями ДУ При этом приращение вектора искомых переменных
ДX=o= ^ДУ
(24)
а вектор приращения уклонений ДV размерности (п-т) в качестве компонентов может содержать как вариации уклонения от тривиальных плоскостей (4), (5), так и от поверхностей, описываемых нелинейными уравнениями (3), аппроксимированные линейной формой (12):
ДV=co1{Дz1,Дvj',Дvk"; Цке Б}.
Пусть в точке Х(1) нарушено Фр (X). Попытаемся в точке X(t+1) уменьшить невязку Фр (X) , не нарушая при этом эффективных в X(t) ограничений. В соответствии с (12),(13) при достаточно малых вариациях ДX
ДФр ^ (X«)^ = арДX (25)
Не теряя общности рассуждений, предположим, что Фр (X(t))< 0, и нам требуется при переходе от X(t)
к
х(‘-
получить
Дфр > 0,
(26)
т.е. сократить величину невязки. Подставляя (29) в (30) и учитывая условие (26), получим:
ДФр = UpДV > 0 (27)
где
ир =(ир1,...,ирг ) = арН
(28)
Отметим специфику случая сокращения невязки равенств: в (28) следует брать ар если Фр)<0, и -ар - если Фр(X(t))>0.
= Аб
и
уз
Таким образом, задача сокращения невязки Фр (X) сведена к выбору допустимого вектора уклонений ДУ , обеспечивающего выполнение условий (27). Аналогичные рассуждения можно провести для этапа улучшения функции цели, где вместо Фр (X) в формулах (25)-(28) будет фигурировать ф0 (X).
Для выполнения условия (27) компоненты ДV=co1(Дvj, ]=1,..,(п-т)) необходимо выбирать по следующим правилам:
Av. =
J
0, если j е 12 ( базисное нелинейное равенство из (3))
0, если jе13(14) 8 ир <0 при Vj(Vj) = 0 (базисноетривиальное ограничение - грань гиперкуба (4),(5))
Как правило, методы, использующие квадратичную аппроксимацию нелинейностей ф. (X), требуют применения матрицы вторых производных (матрицы Гессе) [10], что непомерно увеличивает объем вычислений на каждой итерации.
Поскольку нам необходимо значение второй производной только по направлению о, эту трудность можно обойти, используя приближенное численное дифференцирование вдоль о вместо аналитического расчета Р; через матрицы Гессе.
(29)
u ,если JeT (I.) при v. (vj )>0 ( внутренняя точка гиперкуба (4),(5))
Итак,
Определение длины шага
Пусть в соответствии с (29) вычислен вектор AV. Тогда, в соответствии с (24), можно определить AX и допустимое направление о:
s = col (о4, о2,..., on ) = AX = H*AV (30)
Длина шага (X) вдоль направления о ограничена следующими условиями:
X = min {X', X'', X"'}, (31)
где:
X' = min {Xj, Xj > 0, j £ } - длина шага до j-ой огра-
ничивающей гиперплоскости тривиальных ограничений (4),(5)), не входящей в базис;
X'' = min {X'' Xi > 0,i e } длина шага до возникновения невязок эффективных в X4 нелинейных ба-
dX
d2фi
"dX2"
где
dфi
dX
= Уф(ч -o(t):
фі (X(t) + Хо)-фі (X(t))
X=X(t)
= o(t)V>(t)o(t)‘
X=X(t)
dфi dфi
dX x=x(t+i) dX X=X(t)
X
-фі X(t) + Xo)
x=X(t+1)
X
(34)
(35)
(36)
X — достаточно малый «пробный» шаг вдоль о;
X (‘+1) = X1 + Хо - пробная точка вдоль о :
У2ф(1) — матрица Гессе в точке X(t).
Таким образом, предлагаемый метод сочетает эконо-зисных ограничений (3), превышающих инженерную мичность методов, использующих только первую произ-
точность;
X'" — длина шага до нарушенного ограничения
фр(x(t)) или до максимума ф0 (X(t)) вдоль направления ется по формулам:
водную, с «дальнодействием» квадратичных методов.
В соответствии с (32), длина шага вдоль о вычисля-
о.
Значения минимального шага X до тривиальных границ вычисляется по формуле
*j =
где:
(j -х(‘>)
Х"=
(t), если о ■' >
(t), если 0((t) < 0
(32)
(37)
(ч(і)) - j-ая координата X« ( ДК(4)).
Для вычисления X'', X''' будем использовать квадратичную аппроксимацию:
Дфі (X, о, Х)«аіХ + РіХ2, (33)
где
= dфl
X'" =
dX
X=X(t)
в = і.
Pi 2 dX2
X=X<4
+~, если p = 0 8 ß0 > 0 - ф-ия цели не ограничена сверху! а
---—, если p = 0, а0 > 0,ß0 < 0 - оптимум!
2ß0
а p 2
---—, если pФ 0, => ар < 4ßp9p несовместимость!(тіп несовместимости)
2ßp
9„(X(t))
—p--, если pф 0,ap ф 0,ßp « 0,
(38)
-ap -^ - 4ßp9p(X(t))
2ßp
, если pФ0 и ((9p(X(t))>0,ap >0,ßp <0) или (9p(X(t))<0,ap <0,ßp >0))
-ap + \/ap -4ßprnp(X(t)) p y p------------------, если pФ0 и ((9p(X(t))>0,ap <0,ßp >0) или
2ßp
(9p(X(t))<0,ap >0,ßp <0))
+
p
3
В формуле (37) не фигурирует а;, поскольку для базисных ограничений в точке X(t)
а. = ^ф;0 = 0.
аХ
Алгоритм Этапа II может быть представлен следующими шагами:
Шаг 1. В текущей точке X® вычисление ф;^®), 1еБ(4). Если >£;, - выполнение устранения не-
вязки базисного нелинейного ограничения-равенства, иначе - переход к Шагу 2.
Шаг 2. Нахождение первого неудовлетворенного в текущей точке небазисного |фp(X(t))| >ер, р£ Б®. Если таковых нет, то переход к Шагу 6, иначе - Шаг 3.
Шаг 3. Вычисление Vфp(X) Ех® по формуле (13) и ^ по формуле (28). Вычисление о® по формуле (29),
(30).
Шаг 4. Вычисление X по формулам (31)-(38). Определение небазисной плоскости - г, шаг до которой минимален ^г, V/ или V").
Шаг 5. Вычисление X(t+1) по формуле (24). Вычисление иг по формуле (28) Введение в базис плоскости zг ^г'^г") вместо одной из V/ (V/' ), по которым уклонялись, - выполнение шага жорданова преобразования по формулам (31),(32).Переход к Шагу 1.
Шаг 6. Вычисление Vф0(X(t)) по формуле (26) и ^ по формуле (28). Если ||Vф0(X(t)) || <е0, расчет окончен, иначе - вычисление о® по формуле (29),(30), переход к Шагу 1.
Расчетный пример
Рассмотрим пример расчета элементарной системы пароснабжения:
Рис. 3.
В узел 1 подается пар в количестве Ці тонн/час, из узла 2 потребитель забирает О тонн/час, из узла 3
- оставшиеся 0’ = 04 + 05 + 06 = 01 - 02 тонн/час распределяются на три ветви, ведущие к висящим узлам
4, 5, 6. Необходимо определить минимальный диаметр (максимальное гидросопротивление) трубопровода потока х3, при котором перепад давления на этом трубопроводе Др3 будет в диапазоне
ДР3тіп < Дрз < ДР3тах
а также величины расходов 04, 05, 06 на ветвях, ведущих к висящим узлам.
В первую очередь решим задачу нахождения максимального гидросопротивления потока х3.
Запишем линейные уравнения материального баланса для узлов 1 и 2 (аналог I закона Кирхгофа для электрических цепей) [4]:
У1= -Х1 - Х2 + = 0, (39)
У2= Х2 - Х3 - = 0.
Здесь множество индексов линейных уравнений 11={1,2}.
Далее, для замкнутого контура I (см. рис. 1), можем записать уравнение второго закона Кирхгофа: сумма падений давления в замкнутом контуре равна нулю. В отличие от электрической цепи связь между потоком в ветви X] и падением давления на ней Др; носит нелинейный (квадратичный) характер [4]:
Др] =Р)*Х]*|х]|, )={1,2,3} (40)
Обозначая неизвестное гидросопротивление третьего потока р3 как х4, а Др3 как Х5, запишем два нелинейных уравнения, из которых первое - это определение перепада давления на трубопроводе по формуле (31), а второе - П-ой закон Кирхгофа для контура I:
Zl= +Х4*Х3*|Х3 I - Х5 =0
Z2= -р1*Х1* 1x1 | + р2*Х2* I Х2 I + Х5 = 0, (41)
где:
р1, р2 - известные (заданные) гидросопротивления трубопроводов потоков Х1,Х2,
х1, х2, х3 - искомые значения потоков, х4 - искомое гидросопротивление,
Х5 - искомый перепад давления в потоке Х3.
Здесь множество индексов нелинейных уравнений І2 = {1,2}.
Двусторонние тривиальные ограничения (4), (5) для рассматриваемого примера имеют вид:
у'1= х1 -X1min > 0
^ Х2 -X2min > 0 (42)
v'з= Х3 -X3min > 0 у'4= х4 > 0
у'5= х5 -ДР3тЬ > 0
у"1^-х1 +X1max > 0
-Х2 +X2max > 0 (43)
v"3= -Х3 +X3max > 0 -Х4 +X4max > 0 v"5= -х5 +ДР3тах > 0
Здесь множество индексов тривиальных ограничений ]={1,2,3,4,5}.
Критерий оптимальности в данном случае линеен:
ф0(X)= х4 ^тах (44)
Пусть при расчете гидравлической цепи рис. 1 01=1,0, 02=0,6 , р1=10, р2=1. Тривиальные ограничения (4),(5) выглядят так:
Е
0< хі < 0.8 0< х2 < 0.8 -0.3< х3 <+ 0.3 0 < х4 < 100 -0,5< Х5 < +0,5
Расчет выполняется с точностью е=0,01.
Этап I - введение в базис линейных ограничений-равенств (2).
В соответствии с Шагом 0 Этапа I, сформируем в виде таблицы Штиффеля-Зуховицкого [2] начальную ‘расширенную’ обратную матрицу G/, вектор-строки a1, a2 и рассчитаем вектор-строку ^' по формулам (19),(20):
Таблица 1
(4'),(5')
Начальные значения уклонений и функции цели для этапа 2: Zl=-0.5: z2=-5.31•, фо=0. Промежуточные результаты работы второго этапа представлены в табл. 5:
Таблица 5
№ X Х1 Х2 Х3 Х4 Х5 Zl Z2 Ф0
1 0,45 0,25 0,75 0,15 0 -0,5 -0,5 -0,5625 0
2 0,5625 0,25 0,75 0,15 0 0,0625 -0,0625 0 0
3 0,0667 0,2411 0,7588 0,1588 0 -0,0042 -0,0042 -0,0099 0
4 18,4748 0,2801 0,7198 0,1198 18,4741 0,2587 -0,0064 -0,0082 18,4741
5 20,7797 0,2986 0,7013 0,1013 39,2539 0,3963 -0,0073 -0,0034 39,2539
6 12,1821 0,3058 0,6941 0,0941 51,4360 0,4568 -0,0009 -0,0032 51,4360
7 6,4005 0,3098 0,6901 0,0901 57,8366 0,4791 -0,0090 -0,0045 57,8366
8 4,7344 0,3115 0,6884 0,0884 62,5710 0,4944 -0,0055 -0,0025 62,5710
9 0,2184 0,3116 0,6883 0,0883 62,7895 0,5 -0,0098 -0,0026 62,7895
Базис Хі V!' V2' V/ V/ V5' 1 я1= (-1, -1, 0, 0, 0)
Х1= 1 0 0 0 0 0 я2=( 0, 1, -1, 0, 0)
Х2 = 0 1 0 0 0 0 -0.3 0 -0.5
Х3 = 0 0 1 0 0
Х4 = 0 0 0 1 0
Х5 = 0 0 0 0 1
^'= -1 -1 0 0 0 1.0
Ищем разрешающий элемент - наибольший по модулю коэффициент в вектор-строке (и11=-1, в табл. 1 затенен). Исключаем из базиса V! и вводим в базис у1, выполняя шаг жорданова преобразования по формулам (21),(22).
Получаем новую расширенную обратную матрицу (табл. 2). Для a2 вычисляем находим разрешающий элемент (и22). Исключаем из базиса V/ и вводим в базис у2, выполняя шаг жорданова преобразования. Получаем окончательную обратную матрицу Этапа I (табл. 3).
Таблица 2
Базис Хі У1 V2 ' V/ V/ V5' 1
Х1 = -1 -1 0 0 0 1.0
Х2= 0 1 0 0 0 0
Х3= 0 0 1 0 0 -0.3 0 -0.5
Х4= 0 0 0 1 0
Х5= 0 0 0 0 1
и2'= 0 1 -1 0 0 -0.3
Таблица 3
У1 У2 V/ V/ V5' 1
-1 -1 -1 0 0 0.7
0 1 1 0 0 0.3
0 0 1 0 0 -0.3
0 0 0 1 0 0
0 0 0 0 1 -0.5
Начальная матрица (Н-затенена) и начальный вектор-столбец искомых переменных (Х(0) - светлый столбец) Этапа II изображены в табл. 4:
Таблица 4
Базис Хі V/ V/ V5 ' 1
Х1 = -1 0 0 0.7
Х2 = 1 0 0 0.3
Х3 = 1 0 0 -0.3
Х4 = 0 1 0 0
Х5 = 0 0 1 -0.5
Итеративный алгоритм завершен, так как достигнуто граничное значение для переменной Х5 (Х5 < 0,5), что означает, что мы нашли максимум функции цели при заданных ограничениях.
Таким образом, искомое максимальное гидросопротивления потока х3 равно 62, 79.
Теперь перейдем к задаче расчета расходов О4, О и Об.
Пусть гидросопротивления для соответсвую-щих потоков х4, Х5 и хб равны р4 = 10, р5 = 20, рб = =30.
Для ветви потока х3 выберем гидросопротивление р3 = 57,84 (соответствующий перепад давления Др3 тогда равен 0,48 (согласно табл. 5)). На основании этих величин, рассчитаем значение расхода О3 потока х3 из формулы:
Дк = ^ ■ О2 (45)
Оз2 = Др3 /р3 = 0,48/57,84 = 0,008;
О3 = 0, 09 тонн/час.
Теперь вычислим проводимости для каждой из ветвей:
Уп =1/рп (п=1+к) (46)
У1 =1/р1 =1/10 = 0,1
У2 =1/р2 =1/1 = 1 у3 =1/р3 =1/57,84 = 0,017 У4 =1/р4 =1/10 = 0,1 У5 =1/р5 =1/20 = 0,05 уб =1/рб =1/30 = 0,033
Определим суммарную проводимость для всей сети:
к
Y= £ Уп = У1 + У2 + У3 + У4 + У5 + Уб = 0,
п=1
1+ 1 + 0,017 +0,1 + 0, 05 + 0, 033 = 1,3
Определим суммарную проводимость для ветвей с висящими узлами:
к
^= £ У, = У4 + У5 + Уб = 0,1 + 0, 05 + 0, 033 = 0, 183 1=1
В качестве эквивалентной веуви выберем в£твь потока х3. Определим Р0 = £ У, ■ (Р,+Др,)/£ Уп
1=1 п=1
граничное давление для эквивалентной ветви (Оз(Р0) = 0). Пусть давление Р3 ветви потока х3 фиксировано и равно 1,5 МПа. Тогда, с учетом величины перепада давления, граничное давление для эквивалентной ветви будет равно:
Р0 = £у, (Р + Др,)/£Уп =
1=1 п=1
=0,183 * (1,5 + 0,48) / 1,3 = 0,29
Определим расходы на ветвях с висящими узлами при О3 = 0:
О,,0 = У, ((Р3 +ДР3)-Р0) (47)
О40 = У 4 ■ (Р0 - (Р3 +ДР3)) = 0,1 (1,5+0,48 - 0,29) = 0,17 О50 = у5 ■ (Р0 - (Р3 + Др3)) = 0,05 ■ (1,5 + 0,48 - 0,29) = 0,08 Об0 = уб ■ (Р0 - (Р3 + Др3)) = 0,033 ■ (1,5 + 0,48 - 0,29) = 0,06
На следующем шаге определим доли расхода О3, распределяемых в 1 - ю ветвь:
* = У^Л’ (48)
q4 = у4Л’ = 0,1/0,183 = 0,55 *5 = У5/У’ = 0,05/0,183 = 0,27 *б = Уб/У’ = 0,033/0,183 = 0,18
Теперь рассчитаем эквивалентное сопротивление эквивалентной ветви:
р3,0 = р3 + 1Л’ = 57,84 +1/0,183 = 63,3
И дальше определим расходы в ветвях, содержащих висящие узлы:
О, = О3*1 + О,,0 (49)
О4 = О3*4 + О4,0 = 0,090,55 + 0,17 = 0,22
Об = О3^5 + Об,0 = 0,090,27 + 0,08 = 0,1
О6 = О3* + О6,0 = 0,090,18 + 0,06 = 0,08
Следовательно, искомые расходы О4, Об, О6 в вет-
вях найдены и равны соответственно 0,22; 0,1 и 0,08 тонн/час.
Выводы
В процессе подробного изучения двух современных методов расчета гидравлических сетей, а именно метода топологической свертки и алгоритма решения задач нелинейного программирования с использованием базиса сокращенной размерности - были выделены основные преимущества и сильные стороны этих двух методов.
Важным преимуществом метода топологической свертки является процедура исключения из расчета
висящего узла графа системы. Благодаря этой процедуре, представляется возможным решение проблемы висящих узлов, существующих в больших гидравлических сетях [2].
В то же время, главнейшей особенностью усовершенствованного алгоритма решения задач нелинейного программирования с использованием базиса сокращенной размерности является существенная оптимизация расчета гидравлической сети.
На основании полученных выводов, была установлена возможность эффективного совместного применения вышеуказанных методов для расчета гидравлических сетей с висящими узлами. Идея состоит в том, чтобы использовать процедуру исключения висящих узлов метода топологической свертки и расчетнооптимизационный механизм усовершенствованного алгоритма решения задач нелинейного программирования с использованием базиса сокращенной размерности.
Как показали практические исследования, данный прием особо эффективен для гидравлических сетей сложной структуры с большим количеством висящих узлов.
Литература
1. Щербашин Ю.Д. Метод решения задач нелинейного программирования с использованием базиса переменной размерности. // Проблемы управления и информатики. -2006.-№°3.-С.27-36.
2. Винничук С. Д. Метод линейной свертки для расчета распределительных сетей. // Моделирование и диагностика сложных процессов и систем: Сб. науч. тр. - Киев: ИПМЭ НАН Украины, 1997. - с. 71-79.
3. Гасс С. Линейное программирование. -М.Физматгиз, 1961,-303с
4. Меренков А.П., Хасилев В.Я. Теория гидравлических цепей. -М.: Наука, 1985.- 278с.
5. Дерябин А.Н., Кондратьев А..В., Файзуллин Р.Т., Федоров С.В. Квазистационарный расчет больших электрических цепей с нелинейными элементами. // Труды Международной конференции “Математические модели и методы их исследования”.-: Красноярск - 2001.- Т.1, С.223-227.
6. Щербашин Ю.Д. Об одном методе улучшения непрерывного производственного процесса. // Автоматика и телемеханика. -1970.-№9.-С.119-125.
7. Левин М.М., Волковицкая П.И., Лаптин Ю.П.,Журбенн-ко Н.Г. Использование средств оптимизации в системе автоматизированного проектирования энергетических котлоагрегатов КРОКУС. //Энергетика и электрифика-ция.-2003.-№7.-С41-51.
8. Константинова Е.С., Курдюк Е.В., Щербашин Ю.Д. Оптимизация технологического режима цеха синтеза аммиака. // Сб.”Моделирование и оптимизация каталитических процессов”. -М.:Наука, 1965, -С.109-127.
9. Griffith R. E., Stewart R. A. A non-linear programming technique for the optimization of continuous processing systems.
— Management Sci., 1961, № 7, p.379-392.
E
