Применение динамического программирования при непрерывных переменных

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

Рассмотрим, как решается задача о распределении ресурсов (7.1.1), (7.1.2) при непрерывных переменных. Определим последовательность функций состояния

при условии

.

Основное рекуррентное соотношение имеет вид:

(7.1.19)

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

Рассмотрим задачу определения при фиксированном . Обозначим

.            (7.1.20)

Тогда

.

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

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

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

Заметим, что задача вычисления значительно упрощается при условии, что все выпуклы или вогнуты. Можно показать, что из выпуклости (или вогнутости) всех функций следует выпуклость (вогнутость) функции . Рассмотрим такие случаи.

Случай 1.Все функции выпуклы.Поскольку выпуклая

функция, то и выпуклая относительно переменной . Нас интересует при фиксированном . Как известно, максимум выпуклой функции на замкнутом ограниченном множестве достигается в одной из крайних точек. В данном случае такими точками являются: 1) = 0; 2) = .

Итак, в случае выпуклости всех , для нахождения  достаточно проверить лишь крайние точки = 0 и = .

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

Особенности метода ДП.

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

Характеристика многошаговых распределительных задач.

 

В распределительных задачах с большим числом различных результатов производственной деятельности (i=1,n)  и видов ресурсов (j=1,m) общее решение задачи оптимизации может быть при определенных условиях заменено совокупностью последовательно решаемых менее сложных частных задач оптимизации, например, по каждому из отдельных видов производственной деятельности. При этом важными понятиями ДП являются: последовательность шагов оптимизации, состояние системы распределения ресурсов и варианты решения (области изменения оптимизируемых переменных)

Принципы инвариантного погружения и  оптимальности.

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

z =c₁x₁+c₂x₂+…+c_nx_n max,

a₁₁x₁+a₁₂x₂+…+a_1nx_nb₁,

a₂₁x₁+a₂₂x₂+…+a_2nx_nb₂,

…………………….

a_m1x₁+a_m2x₂+…+a_mnx_nb_m,

x₁,…x_n0,

где z-целевая функция, подлежащая максимизации;

x_i-оптимизируемые переменные;

i=1,n-номер оптимизируемой переменной;

с_i-доход от реализации i-го вида производственной деятельности;

j=1,m-номер ограничений на значения переменных;

а_ij-коэффициенты уравнений-ограничений;

b_j-величина j-го ресурса (правая часть ограничений).

Здесь каждый вид производственной деятельности iможет рассматриваться как отдельныйшагоптимизации; множество возможных значений переменных (допустимая область решений) x_iкакварианты решений, а количество каждогоj-го вида ресурса (B_i1,…, B_ij,… B_im), 0B_ijb_j, доступного для распределения на текущем и предыдущих (либо текущем и последующих) шагах (поi-м типам деятельности) каксостояние модели. Тогда оптимальное значение целевой функцииzдля шагов i,i+1,…,n  при заданных состояниях {B_ij} может быть записано в виде следующей рекуррентной функции Белмана (алгоритма прямой прогонки):

f_i(B_i1,…, B_{i m}) = max  {c_i x_i + f_i-1(B_i1-a_i1x_i ,…,B_{i m}-a_{i m}x_i )}, i =1,n;j =1,m, (1)

                                           0a_ijx_iB_ij

с начальными условиями  f₀(B₀₁,…,B_{0 m})=0.

Оптимальное значение целевой функции в обратном времени для шагов n,…, i, i-1,…,1 при заданных состояниях {B_ij} может быть записано в виде следующего алгоритма обратной прогонки:

f_n(B_n1,…,B_nm)=max{c_nx_n},                                                                                                0a_njx_nB_nj

f_i(B_i1,…,B_im)=max {c_ix_i+f_i+1(B_i1-a_i1x_i,…,B_{i m}-a_{i m}x_i)},i=1,n;j=1,m,  (2 )

                           0a_ijx_iB_ij

где 0B_ijb_j.

Разница между прямым и обратным способами решения задачи заключается в определении состояния модели. В прямой модели B^/_{i j}- количество ресурса j-го типа, распределяемого от первого шага до i-го, а для обратной модели B_ij- количество ресурса, распределяемого на всех шагах от i-го до n-го.

Решение задачи (1) основывается на двух основополагающих принципах.

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

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

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

Методика реализации принципа оптимальности. Примеры.

 

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

Пример№1.Оптимизация производства услуг в сети спутниковой связи.

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

maxz= 3x_a+ 2x_б ,

при ограничениях

                                                                x_а+  2 x_б       6 ,

                                              2 x_а+   x_б 8 ,

                                              - x_а  +  x_б       1,

                                                         x_б 2,

                                                [x_а ] ,[x_б ] 0 .

Здесь шаги оптимизации определяются  порядком оптимизации различных видов услуг связи: на первом шаге i =1 оптимизируется число услуг типа а, на втором шаге i =2 оптимизируется число услуг типаб.

 Вектор состояния (B_i1,B_i2) определяется двумя видами ресурсов: числом телеграфных и телефонных каналов, подлежащих распределению на i-м шаге, т.е.0B_i16 и 0B_i28. Варианты решений определяются допустимой областью определения переменных x_i, подлежащих оптимизации. Далее рассмотрим реализацию алгоритма ДП  для случая обратного алгоритма (1).

Критерий оптимальности для первого шага (оптимизация производства услуги ) принимает вид:

                                            f₂(B₂₁,B₂₂)=max{2x_б},

02x_бB₂₁

0x_бB₂₂

Далее определяем допустимые границы изменения состояния на первом шаге:

Процесс поиска условно оптимальных решений по величине производства услуги x_ботобразим в форме таблицы 1.

   Таблица 1

						Оптимальное решение
	Xб=0	Xб=1	Xб=2	Xб=3	Xб=4
(2,0)	0	-	-	-	-	0	0
(3,2)	0	2	-	-	-	2	1
(4,4)	0	2	4	-	-	4	2
(5,6)	0	2	4	-	-	4	2
(6,8)	0	2	4	6	-	6	3

 

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

Целевая функция на втором шаге имеет следующий вид:

 

 

Решение задачи второго этапа представлено таблицей 2.

 

 

Таблица 2

 

						Оптимальное решение
	Xa=0	Xa=1	Xa=2	Xa=3	Xa=4
(6,8)	0+6=6	3+4=7	6+4=10	9+2=11	12+0	12	4

 

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

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

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

                                            f₂(B₂₁,B₂₂)=max{2x_б}.

02x_бB₂₁

0x_бB₂₂

Так как из ограничений следует, что x_бmin {B₂₁/2,B₂₂}, а  f₂(x_б|B₂₁,B₂₂)= 2x_б, то, подставляя первое во второе и переходя к максимуму целевой функции, получим

f₂(B₂₁,B₂₂) =max { f₂(x_б|B₂₁,B₂₂)=2x_б=

                                          x_б

    =2min{B₂₁/2,B₂₂}, откудаx^*_б=min{ B₂₁/2,B₂₂}.

Т.е. оптимальное значение числа услуг типа б, получаемое на шаге 1 , равно минимальному из двух  видов ресурсов (числу телеграфных B₂₁/2 или телефонных каналов B₂₂), распределяемых на первом шаге.

 Далее для шага 2 имеем:

f₁(B₁₁,B₁₂) = max  {3x_a+f₂(B₁₁-x_a,B₁₂-2x_a)}=

= max  {3x_a+2min{(B₁₁-x_a)/2,B₁₂-2x_a}},

                       0x_aB₁₁;                                                                                                                                                                                    0                              02x_aB₁₂,

где B₁₁=6 и B₁₂=8 для первого шага оптимизации.

Подставляя значение ресурсов ТГ и ТФ каналов в ограничения, получим обобщенное ограничение в виде x_amin(B₁₁,B₁₂/2)=4.                                                                              Учитывая пропорциональную зависимость значения целевой функции от значения x_aпоследнее неравенство превращается в равенство, а соответствующее решению  x_a^*= 4 оптимальное (максимальное) ее значение  f₁^*(4) = 12.

Второй этап. Безусловное решение на первом шаге оптимизации числа предоставляемых услуг x_бнеобходимо проводить для следующего состояния по ресурсам: B₂₁= B₁₁- x_a^*=6-4=2;    B₂₂= B₁₂- 2x_a^*=8-8=0. Откуда из условного решения следует безусловное:

x_б^*= min { B₂₁/2, B₂₂}= min { 1, 0 }= 0 .

Таким образом, пошаговое решение задачи линейного программирования методом динамического программирования обеспечило в целом оптимальную стратегию производства услуг в одну единицу времени  x_a= 4  , x_б= 0 при максимальном для случая целочисленного решения значении дохода сети спутниковой связи z = 12 у.е./ед.вр.

 

Пример№2. Задача оптимального распределения ресурсов резервирования в радиорелейной линии связи.

 

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

где P_i–вероятность безотказной работы i-го интервала.

Для повышения надежности данной последовательной системы используется резервирование станций на каждом из отдельных интервалов РРЛ, вероятность безотказной работы которых определяется выражением:

где q_i  и p_i–вероятность отказа и вероятность безотказной работы элемента  на i-м интервале соответственно; x_i- число резервных станций на i-м интервале; 1+x_i–общее число (одна рабочая и x_iрезервная) станций в i-м интервале.

Пусть также введены ограничения на число резервных станций x_i2;i=1,n. При этом суммарная стоимость резервных элементов с учетом известных ограничений на число станций в подразделении, развертывающем РРЛ не может превысить величины  С= 1200у.ед.

Остальные исходные данные, содержащие сведения о надежности p_i(x_i) и стоимости резервных средств c_i(x_i) i-го интервала, даны в таблице 1.

xi	i=1		i=2		i=3
	p1(xi)	C1(xi)	p2(xi)	C2(xi)	p3(xi)	C3(xi)
0 1 2	0,70 0,91 0,973	100 200 300	0,60 0,84 0,936	200 400 600	0,50 0,75 0,875	300 600 900

 

В таблице 1 значения вероятностей безотказной работы i-х интервалов при использовании в них x_iрезервных станций определяются из выражения

p_i(x_i) = 1-(1-p_i)^1+x_i,

где p_i-вероятность безотказной работы интервала при отсутствии резерва.

Необходимо определить оптимальное количество резервных станций на каждом интервале x_i^опт, обеспечивающих максимальную надежность РРЛ, т.е.

P_л()=

при ограничениях

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

Рекуррентное соотношение для функции Беллмана в данном случае может быть записано в виде:

F_n(s_n) =

F_i (s_i)=

i=1,2,…,n-1.

Из таблицы №1 найдем границы изменения состояния s_iна каждом шаге: s_3min=c₃(x₃=0)=300у.е.,

s_3max=,s_2min=

s_2max,s₁=

На первом этапеобратного алгоритма рассматривают шаг 3. Найденные на этом и последующих этапах значения критерия и условно-оптимальные решения для всех допустимых значений s_iи оптимизируемых переменных x_i=0,1,2 представлены в таблицах №2, 3, 4.

 

Шаг№ 3                                                                          Таблица№2

S3	P3(X3)			Условное оптимальное решение
	X3=0	X3=1	X3=2	F3(s3)	X3^опт
	P=0,5 C3=300	P=0,75 C3=600	P=0,875 C3=900
300 400 500 600 700 800 900	0,5 0,5 0,5 0,5 0,5 0,5 0,5	- - - 0,75 0,75 0,75 0,75	- - - - - - 0,875	0,5 0,5 0,5 0,75 0,75 0,75 0,875	0 0 0 1 1 1 2

 

Шаг№2                                                                     Таблица№3

S2	P2(X2)=p2(x2) F3(s2-c2(x2))			Условное оптимальное решение
	X2=0	X2=1	X2=2	F2(s2)	X2^опт
	P=0,6 C2=200	P=0,84 C2=400	P=0,936 C2=600
500 600 700 800 900 1000 1100	0,3 0,3 0,3 0,45 0,45 0,45 0,525	- - 0,42 0,42 0,42 0,63 0,63	- - - - 0,468 0,468 0,468	0,3 0,3 0,42 0,42 0,468 0,63 0,63	0 0 1 1 2 1 1

 

Здесь значения функции Беллмана на втором шаге определяются с учетом ее значения на предыдущем шаге согласно выражению

F_i(c_i)= p_i(x_i)F_i+1(-c_i).

 

Шаг3                                                                    Таблица№4

S1=	P1(X1)=p1(x1)F2(c2)			Безусловное оптимальное решение
	X1=0	X1=1	X1=2	F1(s1)	X1^опт
	P=0,7 C1=100	P=0,91 C1=200	P=0,973 C1=300
1200	0,441	0,5733	0,4554	0,5733	1

 

Безусловное оптимальное решение на первом этапе получено лишь для шага 3 (числа резервных станций на первом интервале x₁^опт=1), поэтому окончательная стратегия относительно необходимых резервных станций в каждом из трех интервалов может быть получена лишьна втором этапе– этапе анализа результатов пошаговой оптимизации.

Так как  x₁^опт=1, то величина ресурса, подлежащая распределению на первом и втором шагах будет равно S₂=у.е. Тогда из таблицы№3 для S₂=1000 получим оптимальное число резервных станций во втором интервале:  x₂^опт=1. Следовательно, состояние S₃=. Наконец, из таблицы №4 для S₃=600 имеем оптимальное число резервных станций на третьем интервале: x₃^опт=1.

Окончательно, оптимальная стратегия распределения резерва имеет вид: , т.е. на каждом радиорелейном интервале должно находиться по одной резервной радиорелейной станции. Вероятность безотказной работы РРЛ при этом составляет P_РРЛ= 0,5733; Суммарные затраты и затраты по интервалам составляют ; С₁=200; С₂=400; С₃=600у.е.. Интервальные вероятности безотказной работы составляют p₁(x₁=1)=0,91; p₂(x₂=1)=0,84; p₃(x₃=1)=0,75.

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

(17_05)

Еще раз рассмотрим постановку задачи динамического программирования рассмотрим на примере инвестирования, связанного с распределением средств между предприятиями. В результате управления инвестициями система последовательно переводится из начального состояния S₀В конечноеS_n. Предположим, что управление можно разбить на n шагов и решение принимается последовательно на каждом шаге, а управление представляет собой совокупность n пошаговых управлений. На каждом шаге необходимо определить два типа переменных - переменную состояния системыS_kпеременную управленияx_k . ПеременнаяS_kопределяет, в каких состояниях может оказаться система на рассматриваемомk-м шаге. В зависимости от состоянияSна этом шаге можно применить некоторые управления, которые характеризуются переменной xk которые удовлетворяют определенным ограничениям и называются допустимыми. Допустим.= (x₁ , x₂ , ..., x_k , ..., x_n )- управление, переводящее систему на состоянияS₀в состояниеS_n , aS_k- есть состояние системы наk-м шаге управления. Тогда последовательность состояний системы можно представить в виде графа, представленного на рис. 2.11.

рис. 2.11

Применение управляющего воздействия x_k на каждом шаге переводит систему в новое состояниеS¹  (S, x_k)и приносит некоторый результатW_k (S, x_k). Для каждого возможного состояния на каждом шаге среди всех возможных управлений выбирается оптимальное управлениеx*_k, такое, чтобы результат, который достигается за шаги сk-го по последнийn-й, оказался бы оптимальным. Числовая характеристика этого результата называется функцией БеллманаF_k (S)и зависит от номера шагаkи состояния системыS.

Задача динамического программирования формулируется следующим образом: требуется определить такое управление x_, переводящее систему из начального состояния S₀в конечное состояниеS_n, при котором целевая функция принимает наибольшее (наименьшее) значениеF(S₀ ,x_) => extr.

Рассмотрим более подробно особенности математической модели динамического программирования:

задача оптимизации формулируется как конечный многошаговый процесс управления;

целевая функция (выигрыш) является аддитивной и равна сумме целевых функций каждого шага:

выбор управления x_kна каждом шаге зависит только от состояния системыkэтому шагуS_k-1, и не влияет на предшествующие шаги (нет обратной связи);

состояние системы S_kпосле каждого шага управления зависит только от предшествующего состояния системыS_k-1и этого управляющего воздействияx_h(отсутствие последействия) и может быть записано в виде уравнения состояния:S_k = f_k(S_k-1 , x_k), k = 1, n;

на каждом шаге управление x_kзависит от конечного числа управляющих переменных, а состояние системы зависитS_k– от конечного числа параметров;

оптимальное управление представляет собой вектор , определяемый последовательностью оптимальных пошаговых управлений:= (x*₁ , x*₂ , ..., x*_k , ..., x*_n ), число которых и определяет количество шагов задачи.

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

Данный пример демонстрирует следующий факт: в многошаговых процессах управление на каждом конкретном шаге надо выбирать с учетом его будущих воздействий на весь процесс. Кроме того, при выборе управления на данном шаге следует учитывать возможные варианты состояния предыдущего шага. Например, при определении количества средств, вкладываемых в предприятие в i-м году, необходимо знать, сколько средств осталось в наличия к атому году и какой доход получен в предыдущем (i- 1)-м году. Таким образом, при выборе шагового управления необходимо учитывать следующие требования:

возможные исходы предыдущего шага S_k-1 ;

влияние управления x_kна все оставшиеся до конца процесса шаги (n-k).

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

Условная оптимизация

На первом этапе решения задачи, называемом условной оптимизацией, определяются функция Беллмана и оптимальные управления для всех возможных состояний на каждом шаге, начиная с последнего в соответствии с алгоритмом обратной прогонки. На последнем, n-м шаге оптимальное управление -х*_nопределяется функцией Беллмана:F(S) = max {W_n (S, x_n )}, в соответствии с которой максимум выбирается из всех возможных значенийx_n, причемx_n € X. Дальнейшие вычисления производятся согласно рекуррентному соотношению, связывающему функцию Беллмана на каждом шаге с этой же функцией, но вычисленной на предыдущем шаге. В общем виде это уравнение имеет видF_n(S) = max {W_n (S, x_n ) + F_k+1(Sⁿ(S, x_k )}x_k € X. Этот максимум (или минимум) определяется по всем возможным дляkиSзначениям переменной управленияX.

Безусловная оптимизация

После того, как функция Беллмана и соответствующие оптимальные управления найдены для всех шагов с n-го по первый, осуществляется второй этап решения задачи, называемый безусловной оптимизацией. Пользуясь тем, что на первом шаге (k = 1) состояние системы известно - это ее начальное состояниеS₀ , можно найти оптимальный результат за всеnшагов и оптимальное управление на первом шагеx₁, которое этот результат доставляет. После применения этого управления система перейдет в другое состояниеS¹(S, x*₁ ), зная которое, можно, пользуясь результатами условной оптимизации, найти оптимальное управление на втором шагеx*₂ , и так далее до последнегоn-го шага. Вычислительную схему динамического программирования можно строить на сетевых моделях, а также по алгоритмам прямой прогонки (от начала) и обратной прогонки (от конца к началу).

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

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

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

 

Рис. 2.12

В задаче имеется ограничение - двигаться по изображенным на схеме маршрутам можно только слева на право, т.е. попав, например, в пункт 7, мы имеем право переместиться только в пункт10и не можем возвратиться обратно в5-й или6-й. Эта особенность транспортной сети дает право отнести каждый из десяти пунктов к одному из поясов. Будем считать, что пункт принадлежитk-му поясу, если из него попасть в конечный пункт ровно заkшагов, т.е. с заездом ровно в (k - 1)-й промежуточный пункт. Таким образом, пункты7, 8 и9 принадлежат к первому поясу,5и6 - ко второму,2, 3и4- к третьему и1- к четвертому. Тогда наk-м шаге будем находить оптимальные маршруты перевозки груза из пунктовk-го пояса до конечного пункта. Оптимизацию будем производить с конца процесса, и потому, дойдя доk-го шага, неизвестно, в каком из пунктовk-го пояса окажется груз, перевозимый из первого пункта.

Введем обозначения: k- номер шага (k = 1, 2,3,4);i- пункт, из которого осуществляются перевозки (i = 1,2,..., 9); j- пункт, в который доставляется груз (j = 2,3,.., 10); С_{i, j}- стоимость перевозки груза из пункта iв пункт j.F_k (i)- минимальные затраты на перевозку груза наk-м шаге решения задачи из пунктаiдо конечного пункта.

Очевидно, что минимум затрат на перевозку груза из пунктов k-го пояса до пункта10будет зависеть от того, в каком пункте этого пояса мы оказались. Номер iпункта, принадлежащегоk-му поясу, будет являться переменной состояния системы наk-м шаге. Поскольку оптимизация осуществляется с конца процесса, то, находясь в некотором пунктеik-го пояса, принимается решение о перемещении груза в один из пунктов (k– 1)-го пояса, а направление дальнейшего движения известно из предыдущих шагов. Номер jпункта (k- 1)-го пояса будет переменной управления на k-м шаге.

Для первого шага управления (k - 1) функция Беллмана представляет собой минимальные затраты на перевозку груза из пунктов1-го пояса в конечный пункт, т.е.F₁(i)=С_{i  10}. Для последующих шагов затраты складываются из двух слагаемых - стоимости перевозки грузаС_{i, j}из пунктаi k-го пояса в пункт j(k- 1)-го пояса и минимально возможных затрат на перевозку из пункта j до конечного пункта, т.е.- F_{k - 1} (i). Таким образом, функциональное уравнение Беллмана будет иметь вид

Минимум затрат достигается на некотором значении j*, которое является оптимальным направлением движения из пунктаiв конечный пункт. На четвертом шаге попадаем на 4-й пояс и состояние системы становится определенным i = 1. ФункцияF₄(1)представляет собой минимально возможные затраты по перемещению груза из1-го пункта в10-й. Оптимальный маршрут определяется в результате анализа всех шагов в обратном порядке, а выбор некоторого управленияjнаk-м шаге приводит к тому, что состояние системы на (k- 1)-м шаге становится определенным.

Решим сформулированную выше задачу, исходные данные которой приведены на рис. 2.12

I этап.Условная оптимизация. 1-й шаг.k = 1F₁(i)=С_{i  10}На первом шаге в пункт10груз может быть доставлен из пунктов7,8или9.

Таблица. 2.18

2-й шаг. k = 2 Функциональное уравнение на втором шаге принимает вид

Все возможные перемещения груза на втором шаге и результаты расчета приведены в табл. 2.19

Таблица 2.19

3-й шаг. k = 3.

Таблица 2.20

4-й шаг. k = 4.

Таблица 2.21

II этап. Безусловная оптимизация

Рис.2.13.

На этапе условной оптимизации получено, что минимальные затраты на перевозку груза из пункта 1в пункт10составляютF₄(1) = 20. Данный результат достигается при движении груза из1-го пункта в3-й. По данным табл. 2.20, из пункта3необходимо двигаться в пункт6, затем - в пункт7(см. табл.2.19) и из него - в конечный пункт (см. табл. 2.18). Таким образом, оптимальный маршрут доставки груза:1=>3=>6=>7=>10. (На рис.2.13 он выделен толщиной стрелок.)

Лекция. .

Формулировка метода конечных элементов в глобальных координатах для двумерной задачи теплопроводности

Вариационные задачи

К решению встречающихся в технике задач прикладной

механики существуют два подхода . При одном из них

используются дифференциальные уравнения, описывающие поведение

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

подход состоит в том, что постулируется вариационный

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

решение минимизирует некоторую величину \chi, которая

определяется как некоторый интеграл от неизвестных величин по

всей области. Интегральную величину \chi, представляющую собой

функцию от неизвестных функций, называют функционалом.

С математической точки зрения оба эти подхода

эквивалентны, и решение, полученное при одном подходе, является

решением при другом подходе. Каждый из этих подходов

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

первый. От одного подхода можно перейти к другому с по-

мощыо математических преобразований, что является

предметом многочисленных книг по вариационным методам .

Различие между этими подходами состоит в способах

получения приближенного решения. Конечно-разностные

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

разностными; метод Ритца и метод конечных

элементов— связаны с приближенной минимизацией функционала.

Пусть физическая (или чисто математическая) постановка

задачи требует минимизации функционала \chi в некоторой

области. Величина \chi определяется в виде интеграла по области V

и части границы S, на которой неизвестны функция {ф} или

ее производные, т. е. она имеет вид

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

части, подобласти, которые будем далее называть элементами,

и пусть функции, которые мы хотим определить, для каждого

элемента записываются в виде

{ф}=[N]{Ф}^e.

Здесь {Ф}^е может содержать узловые значения функции,

соответствующие такому элементу, или некоторые

характеризующие его параметры. Неизвестная функция взята в фигурные

скобки, чтобы показать, что она может быть вектором, a [N]—матрица, определяющая зависимость функции формы от координат.

Для минимизации функционала \chi по всем параметрам {Ф}

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

Если справедливо утверждение, что функционал равен

сумме вкладов отдельных элементов, т. е. что

Х=\sum{}х^e,

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

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

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

минимизирующих функпионал, для всего ансамбля.

Рассматриваем такой подход для плоской задачи теплопроводности.

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

Рассмотрим двумерную задачу теплопроводности через брус квадратного сечения (рис. 1). На верхнем торце бруса поддерживается температура 100оС, на нижнем 50оС, а боковые поверхности идеально изолированы. Требуется найти распределение температуры в брусе, и в частности температуру в точке A (рис. 1). Воспользуемся вариационной формулировкой метода конечных элементов, в которой применяется глобальная система координат Оxy.

Рис. 1. Двумерная задача теплопроводности в брусе квадратного поперечного сечения.

В данном случае определяющим уравнением является уравнение Лапласа.

в D (1)

с граничными условиями Дирихле на части границы:

Т=50, y=0, (2a)

Т=100, y=L (2б)

и условиями Неймана на остальной части границы:

x=0, (3а)

x=L. (3б)

Задача определена на множестве R, состоящем из области D и границы S, т.е. R=D+S.

Уравнения (1) – (3) не будут использоваться непосредственно. Вместо них будет построена эквивалентная вариационная формулировка. С помощью вариационного исчисления можно показать, что решение Т(x, y), удовлетворяющее уравнениям (1) – (3), совпадает с функцией, которая минимизирует функционал.

(4)

Где - функция из допустимого множества пробных функций, заданных в D.

Для этой задачи пробные функции являются допустимыми , если они непрерывны и имеют кусочно-непрерывные первые производные . Кроме того, пробные функции должны удовлетворять главным граничным условиям (2.2). Граничные условия Неймана (3) будут выполняться автоматически для функции, минимизирующей функционал (4), как естественное следствие вариационной формулировки и, следовательно, будут называться естественными граничными условиями.

Разобьем область на конечных элементов. В рассматриваемом примере в качестве конечного элемента выбран треугольник (рис. 2). Общее чисто узлов обозначим n. В двух- и трехмерном случаях нет определенной связи между общим числом элементов и общим числом узлов. Разбиение области и условия непрерывности, накладываемые на пробные функции, позволяют записать функционал (4) в виде

(5)

где - элементный вклад, определяемый равенством

(6)

Рассмотрим типичный элемент ei, показанный на рис. 3. Номера узлов i, j и m должны быть указаны в порядке, соответствующем движению против часовой стрелки.

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

x, yei , (7)

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

С целью определения этих постоянных запишем последовательно уравнение (7) для узлов i, j и m:

=1+2xi+3yi, (8a)

=1+2xi+3yi, (8б)

=1+2xm+3ym, (2.8в)

Рис. 2. Разбиение области на конечных элементов.

Рис. 3. Типичный треугольный элемент e_i .

где ,и- значения Т в узлах i, j и m соответственно. Верхний индекс еi опущен ради упрощения записи. Система уравнений (8) имеет единственное решение для постоянных α1, α 2 и α 3, так как определитель ее матрицы коэффициентов не равен нулю, т.е.

(9)

С использованием тригонометрии легко установить, что этот определитель равен удвоенной площади треугольника, как это и показано в (9). Так как площадь треугольника никогда не равна нулю, т.е. , то решение α 1, α 2 и α3 существует и единственно. Решая (8), получим для α 1, α 2 и α 3 следующие выражения:

(10а)

(10б)

(10в)

где

(11)

остальные постоянные aj,, am, bj, bm, cj, cm могут быть определены путем циклической перестановки индексов. В приведенных выше выражениях для a, b, c и верхний индекс ei вновь опущен, чтобы не усложнять запись формул. Подстановка выражений (10) в (7) дает следующее представление через базисные функции:

(12а)

или

(12б)

где

(13а)

Является матрицей базисных функций, а

(13б)

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

Требуемые производные можно получить, дифференцируя (12а):

(14а)

(14б)

Подстановка (14) в выражение для элементного вклада (6) дает

(15)

Так как подынтегральное выражение в (15) не зависит от x и y и, кроме того,

, (16)

то равенство (15) может быть переписано следующим образом:

. (17)

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

(18)

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

p=1, 2,…,n. (19)

Подстановка (5) в (19) позволяет представить эти уравнения следующим образом:

р=1, 2, …, n. (20)

Очевидно, что при суммировании в (20) ненулевой вклад дают только те элементы, которые содержат узел р. Дифференцирование равенства (17) по позволяет определить вкладэлемента еi в выражении (20). Таким образом, если номера узлов i, j и m элемента еi имеют во множестве номеров узлов системы значения p, q и r соответственно, то дифференцирование (17) поприводит к выражению

(21)

Объединение компонент элементных уравнений, задаваемое равенством (20), называется объединением по узлам, так как процесс объединения должен быть выполнен отдельно для каждого узла системы. Формирование вкладов и их объединение будет проиллюстрировано ниже. На рис. 4 показано разбиение области рассматриваемой задачи на 16 элементов с общим числом узлов, равным 15.

В таб. 1 указаны номера узлов i, j и m для каждого элемента системы. На рис. 4 номера узлов приведены согласно правилу обхода элемента против часовой стрелки. Размеры всех элементов одинаковы, и имеется только два типа элементов, отличающихся друг от друга порядком расположения узлов. Это не является обязательным требованием, но упрощает дальнейшие вычисления.

Рис. 4. Разбиение области на 16 конечных элементов.

Координаты узлов даны в табл. 2. В табл. 3 приведены значения параметров bi, bj, bm, ci, cj и сm, вычисленные по формулам (11). Для иллюстрации рассмотрим элемент 5. Из табл. 1 видно, что его узлы i, j и m имеют номера 7, 4 и 8 соответственно. Подставляя эти данные и соответствующие значения параметров для элемента 5 из табл. 3 в равенство (17), получим для элементного вклада выражение

(22)

Таблица 1

Соотношение между глобальными и локальными номерами узлов

Элемент	Глобальный номер узла			Элемент	Глобальный номер узла
	i	j	m		i	j	m
1 2 3 4 5 6 7 8	4 2 5 3 7 5 8 6	1 5 2 6 4 8 5 9	5 1 6 2 8 4 9 5	9 10 11 12 13 14 15 16	10 8 11 9 13 11 14 12	7 11 8 12 10 14 11 15	11 7 12 8 14 10 15 11

Таблица 2

Координаты узлов

Узел	Координаты		Узел	Координаты		Узел	Координаты
	x	y		x	y		x	y
1 2 3 4 5	0 1 2 0 1	0 0 0 0.5 0.5	6 7 8 9 10	2 0 1 2 0	0.5 1 1 1 1.5	11 12 13 14 15	1 2 0 1 2	1.5 1.5 2 2 2

Таблица 3

Характерные параметры элементов

Элемент	Параметры
	bi	bj	bm	ci	cj	cm
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16	-0,5 0,5 -0,5 0,5 -0,5 0,5 -0,5 0,5 -0,5 0,5 -0,5 0,5 -0,5 0,5 -0,5 0,5	0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0	-0,5 0,5 -0,5 0,5 -0,5 0,5 -0,5 0,5 -0,5 0,5 -0,5 0,5 -0,5 0,5 -0,5 0,5	1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1	1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1	0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

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

i=1, 2,…, 16. (23)

Из выражения (22) видно, что производная не равна нулю только при p=7, 4 или 8. Другими словами, элемент 5 дает вклад только в уравнения системы, включающиеи. Согласно (22), ненулевые производные для элемента 5 определяются выражениями

(24а)

(24б)

(24в)

Уравнения (24) в матричной форме образуют элементное матричное уравнение:

(25)

Три производные относительно ив уравнениях (24) были специально записаны в порядке, соответствующем порядку узлов i, j и m для этого элемента (табл. 1). Таким образом, уравнение (25) – это частный случай (для элемента 5) общего элементного матричного уравнения

(26а)

или

(26б)

где

(27)

есть элементная матрица жесткости элемента e и Тe- элементный узловой вектор, определяемый равенством (13б). Для упрощения записи в дальнейшем верхний индекс ei в уравнениях (26) и (27) опускается.

Так как элементы 1, 3, 5, 7, 9, 11, 13, 15 имеют одинаковые размеры и ориентацию по отношению к системе отсчета Оху, то можно показать, что элементная матрица жесткости k для этих элементов одна и та же при условии, что уравнение (26а) подставлены соответствующие номера узлов. Так, для элемента 9 матричное уравнение имеет вид

(28)

Остальные элементы (2, 4, 6, 8, 10, 12, 14, 16) также имеют общую элементную матрицу жесткости k, которая в этой задаче (что можно проверить) тождественна полученной матрице жесткости для элементов с нечетными номерами. Так, для элемента 10 матричное уравнение имеет вид

(29)

После определения вкладов для всех элементов и объединения их получим матричное уравнение системы.

Для объединения по узлам основным объединяющим уравнением является уравнение (20). Например, для узла 7 соответствующим объединяющим соотношением является

(30а)

Из рис. 4 видно, что для узла 7 дают вклад лишь элементы 5, 9 и 10, поэтому выражение (30а) сводится к

(30б)

Вклады элементов 5, 9 и 10, входящие в правую часть равенства (30б), уже были вычислены; они определяются равенствами (25), (28) и (29) соответственно. Подставляя соответствующие выражения в (30б), получим (31)

и после получения подобных членов

(32)

Уравнение (32) можно также записать в расширенной форме:

(33)

где узловой вектор Т системы имеет вид

(34)

Для каждого узлового параметра системы можно получить уравнение аналогичное уравнению (33). Объединение этих уравнений в одно дает матричное уравнение системы:

(35)

его необходимо подправить для учета граничных условий Дирихле. Введение граничных условий Дирихле приводит к следующему матричному уравнению системы:(36)

Используя любую стандартную процедуру, можно найти что решением уравнения (36) является

(37)

Следовательно, требуемым решением в точке А будет

(38)

Отметим, что здесь не предпринималось каких-либо специальных мер для того, чтобы учесть граничные условия Неймана (3а) и (3б) в узлах 4, 7, 10 и 6, 9, 12 соответственно, поскольку, как ранее установлено, минимизация в этих граничных точках является достаточной для выполнения граничных условий Неймана в качестве естественного следствия.

Переформулируем теперь задачу, которую иллюстрирует рис. 1, с использованием локальной системы координат вместо примененной ранее глобальной системы координат. Вберем квадратичную аппроксимацию вместо линейной пробной функции. Возьмем шестиузловой треугольный элемент (рис 5) с тремя узлами в вершинах и тремя узлами в серединах сторон треугольника, используя в качестве узловых параметров только значения функции в каждом из узлов. Числа 1,2,…,6 обозначают номера узлов такого элемента, занумерованных в соответствии с обходом треугольника против часовой стрелки.

В глобальных координатах XиYполином второго порядка, являющийся пробной функцией для шестиузлового треугольного элемента,

имеет вид

=α+ αx + αy + αx+ αxy + αy, (39)

где α, α, … α- постоянные , в общем случае свои для каждого элемента. Для простоты записи верхний индексв правой части выражения (39) опущен.

Рис 5. Шестиузловой треугольный элемент

Рис. 6. Локальная система координат шестиузлового треугольного элемента.

Если используется локальная система координат ōζη, показанная на рис. 6, то квадратичная пробная функция может быть записана следующим образом:

=α+ αξ + αη + αξ+ αξη + αη, (40)

Здесь также опущены верхние индексы у α .

Далее будет показано, что характерные размеры a, b, c треугольника (рис. 6) играют существенную роль при формулировке задачи в локальной системе координат. Их значения могут быть вычислены по известным глобальным координатам узлов 1, 3, 5 следующим образом. Из тригонометрических соображений получаем

r= [(x₃–x₁)²+ (y₃–y₁)²]^1/2, (41)

cosθ= (x₃–x₁)/r, (42а)

sinθ= (y₃–y₁)/r, , (42б)

где r – расстояние между точками 1 и 3. Характерные размеры a, b, c вычисляются по формулам

a = (x₃ – x₅)* (x₃ – x₁) - (y₅ – y₃)*sin θ ,

b = (x₅ – x₁)* (x₃ – x₁) + (y₃ – y₁)*sin θ , (43)

c = (y₅ – y₃)* (x₃ – x₁) + (x₃ – x₅)*sin θ ,

подстановка формул (42) в (43) дает следующие выражения для характерных размеров:

a = [(x₃ – x₅)*cos θ - (y₅ – y₃)* (y₃ – y₁)]/r ,

b = [(x₅ – x₁)*cos θ + (y₃ – y₁ )* (y₃ – y₁)]/r , (44)

c = [(y₅ – y₃ )*cos θ + (x₃ – x₅)* (y₃ – y₁)]/r ,

которые и являются требуемыми соотношениями.

Для определения постоянных α_iприменим формулу (40) к каждому из узлов 1, 2, …,6 по очереди. После подстановки узловых координат ξ, η, выраженных через величины a, b, c, результирующие уравнения могут быть записаны в матричной форме

(45

Уравнение (45 может быть решено относительно α_iпри условии, что определитель матрицы коэффициентов не равен нулю. Можно показать что определитель равен выражению –c⁴(a + b)⁴/64, которое никогда не обращается в нуль, так как площадь треугольника, равная c*(a + b)/2, никогда не равна нулю. Следовательно, матрица коэффициентов уравнения (45) не сингулярна и имеет обратную.

Обозначая символом А матрицу коэффициентов уравнения (45), можем записать это уравнение в виде

Аа = T(46)

где

=[]=, T=, (47)

Причем Tобозначает узловой вектор элемента. Умножение уравнения (46) на матрицу, обратную А, дает

=АT=B T, (48)

А=B=[b] (49)

Следующим шагом в рассматриваемой здесь формулировки является представление элементного вклада, определяемого равенством (6), через локальные координаты ξ и η. Из (6) видно, что эквивалентные выражения через локальные координаты требуются для производных /x и/y.

Можно показать, что перенос начала координат не изменяет ни узловых значений, ни узловых производных; это значит, что они сохраняют свое численное значение как в глобальной, так и в локальной системе координат. При повороте системы координат узловые значения по-прежнему не меняются, но производные (такие, как /x,/ξ) не являются одинаковыми в двух различных системах координат. Так как пробная функция в равенстве (40) является результатом интерполяции узловых значений,, …[уравнение (45)] и не содержит каких- либо производных, то элементный узловой вектор Tодинаков как в глобальной, так и в локальной системах координат. Следовательно, в этом примере нет необходимости использовать нижние и верхние индексы для обозначения используемой системы отсчета.

Рассмотрим теперь соотношения между координатами x, y и ξ, η некоторой точки. С использованием элементарной тригонометрии эти соотношения, получающиеся в результате поворота на угол , могут быть записаны в матричной форме:

=(50.а)

Обратные соотношения имеют вид

=(50.б)

Квадратная матрица в правой части равенства (50.б) называется матрицей вращения. Ее мы будем обозначать символом R. Обратная матрице вращения матрица Rвходит в равенство (50.а). Можно показать, что матрица, обратная R, равна ее транспонированной матрице это верно для любых ортогональных преобразований).

Если рассматривается как функция ξ, η, т.е.=( ξ, η), то можно записать

=+, (51)

=+

Дифференцирование равенства (50.б) приводит к соотношениям

=cos,=sin, (52)

=-sin,=cos,

а подстановка (52) в (51) дает

= cos- sin, (53)

= sin+ cos

Или в матричных обозначениях

=(54.а)

Обратным равенству (49.а) будет

=(54.б)

Равенства (54) позволяют преобразовать глобальные первые производные в локальные, и наоборот. Необходимо отметить появление в равенстве (54.б) матрицы вращения, а в равенстве (54.а) – обратной ей.

Теперь подынтегральное выражение в равенстве (6) может быть записано через локальные координаты посредством подстановки соотношений (53) или (54.а). Выполнив это, получим

+=+(55)

На основе математического анализа имеем:

=(56)

где и-эквивалентные выражения для функции в двух возможных системах координатxy иξη, а-якобиан преобразования(x, y)/(ξ, η). Из формулы (56) для заменыв выражении (6) имеем

====

=(58)

Здесь использовано выражение (50.а) для производных ,,,. Подстановка (55) и (57) в (6) дает следующее выражение для элементного вклада

=(58)

Представление пробной функции (40) в виде

=(59)

(Здесь используются …и…из табл. 4) позволяет записать производные/,/следующим образом

(60)

Подставляя выражения (60) в (58) получим следующие выражения для элементного вклада:

Введение обозначения

=(62)

Таблица 4

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

1 2 3	0 1 0	0 0 1	4 5 6	2 1 0	0 1 2

Позволяет записать (61) в виде

(63)

В матричной форме равенство (63) имеет вид

(64)

Где определяется выражением (47), а

(65)

Является интегрирующей матрицей. Подстановка (48) в (64) дает формулу для элементного вклада через узловые значения:

(66)

(67)

где (68)

оказывается элементной матрицей жесткости элемента.

Дифференцирование равенства (67) дает элементное матричное уравнение

(69)

Интегрирующая матрица еще должна быть определена явно. Для вычисления элементовнеобходимо [уравнение (62)] вычислить интегралы вида

(70)

Где -целые числа. С учетом рис. 4 этот интеграл может быть записан следующим образом

+(71)

Где первый и второй интегралы в правой части равенства являются вкладами, соответствующими подтреугольникам I и II .

Замена переменной

(72)

Позволяет записать первый интеграл в виде

(73)

Или после интегрирования по

(74)

Формула (73) может быть преобразована к виду

(75)

С использованием замены переменной

(76)

Интеграл в выражении (75) является бета - функцией , которая может быть выражена через гамма – функции и факториалы:

==(77)

Подстановка выражения (77) в (75) дает

=(78)

Аналогично может быть показано, что

=

Подставляя выражение (78) в (70), получим

=(79)

Вычисляя выражение

, (80)

получим искомое значение (62).

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

В таблице 5 выписаны координаты узлов xиy. В таблице 6 дана связь между локальными и глобальными номерами узлов, а также указаны характерные размерыa,bиcдля каждого элемента.

Подстановка величин a,bиcиз таблицы 6 в уравнение (45) приводит к следующей матрице коэффициентов А:

А=(81)

Обратная ей матрица В имеет вид

В=А(82)

Рис. 7. Разбиение области на восемь конечных элементов.

Таблица 5

Координаты узлов

Узел	Координаты		Узел	Координаты			Узел	Координаты
	x	y		x		y		x	y
1 2 3 4 5 6 7 8 9	0 0.5 1 1.5 2 0 0.5 1 1.5	0 0 0 0 0 0.5 0.5 0.5 0.5	10 11 12 13 14 15 16 17 18	2 0 0.5 1 1.5 2 0 0.5 1	0.5 1 1 1 1 1 1.5 1.5 1.5		19 20 21 22 23 24 25	1.5 2 0 0.5 1 1.5 2	1.5 1.5 2 2 2 2 2

Таблица 6

Соотношение между глобальными и локальными номерами узлов и характерные размеры элементов.

Элемент	Номер узла										Размеры
	1	2		3		4		5		6	a	b	C
1 2 3 4 5 6 7 8	11 3 13 5 21 13 23 15	6 8 8 10 16 18 18 20	1 13 3 15 11 23 13 25		7 7 9 9 17 17 19 19		13 1 15 3 23 11 25 13		12 2 14 4 22 12 24 14		1 1 1 1 1 1 1 1	0 0 0 0 0 0 0 0	1 1 1 1 1 1 1 1

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

,=1,2,…,8. (83)

Наконец, подстановка выражения (83) в (80) дает требуемые элементы .

В программе для ЭВМ величины вычисляются непосредственно с использованием формулы (79) и известных значений. Нулевые значенияимогут привести к ошибке, если в программе не позаботиться об этом заранее. Рассмотрим, например, случайи. Из табл. 1 следует, что,,. Подставляя эти значения вместе сив выражение (79) дляи, получим

(84)

Попытка непосредственно вычислить выражение (84) в данном случае привел бы к сообщению об ошибке в программе, так как член [здесь он получается равным (-0)] не определен. Эта трудность возникает когда выражениеимеет минимальное значение, равное -1, т.е. когда. Как видно из формулы (80), выражение (84) не нужно вычислять, если произведениеравно нулю, т.е. либо, либоравно нулю. Таким образом, эта трудность может быть преодолена посредством проверки на равенство нулю произведенияи вычисления членатолько в том случае, когда это произведением не равно нулю. Аналогично величинавычисляется тоько в том случае, когдане равно нулю.

Применяя описанную выше процедуру при вычислении выражения (80), получим матрицу

G=(85)

которая в результате вычислений по формуле (84) принимает вид

G=(86)

Вычисление произведения (68) с использованием (82) и (84) дает следующее выражение для элементной матрицы жесткости:

=(87)

Очевидно, что эта матрица симметрична. Можно доказать, что матрица при заданном порядке локальной нумерации узлов в выражении (87) одинакова для всех элементов. Замена локальных номеров узлов глобальными позволяет объединить матричные уравнения элементов в матричные уравнения системы. Можно проверить, что матричным уравнениям системы после учета граничных условий Дирихле, является уравнение (88). Невыписанные элементы матрицы коэффициентов в этом уравнении равны нулю.

(88)

Решая уравнение (88) каким-либо стандартным способом, получим следующие узловые значения:

,,,,

,,,,,

,,,,, (89)

,,,,,

,,,,,

и, наконец, находим требуемое решение в точке А (рис. 5):

|=87,5°С. (90)