Дворецкий - Основы проектирования химических производств - 2005

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

g	j	(d,ui, ξi) ≤ 0, j J, i I	1	(4.37)
и
				(4.38)
		Bep[ξ Ξ ]≥ ρзад .

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

Нами разработан алгоритм решения задачи (4.35) – (4.38), базирующийся на методе имитационного моделирования [41].

Задача 3. Имеются конструктивные и управляющие переменные. Вектор неопределенных параметров состоит из двух подвекторов ξ1 и ξ2 (ξ = (ξ1, ξ2)). В подвектор ξ1 входят параметры, которые мо-

гут быть только определены на стадии эксплуатации процесса, в подвектор ξ2 – параметры, имеющие неопределенности на этапе эксплуатации те же, что и на этапе проектирования. Пусть при этом ξ1 Ξ1 и

ξ2 Ξ2 .

Эта задача в большей степени соответствует реальным задачам проектирования, поскольку внешние случайные факторы всегда будут иметь место не только на стадии проектирования, но и на стадии эксплуатации производства. Математическая постановка задачи имеет вид:

						Bepξ2 [gj (d,u, ξ1, ξ2) ≤ 0]≥ ρзад,
				(C(d,u, ξ1, ξ2)
min ∫	min M		ξ2				j J	×
d		u
Ξ

×P(ξ1)dξ1 + ∫ minu Mξ2 (C(d,u,ξ1,ξ2)+

Ξ\Ξ

	(gj (d,u, ξ1, ξ2) ≤ 0));
+ A max max(ρзад − Bepξ2		0 P(ξ1)dξ1,
j J*
			:min max [ρзад − Bepξ2 [gj (d,u, ξ1, ξ2) ≤ 0]≤ 0,					(4.39)
						ξ1		, (4.40)
	Ξ = ξ1						Ξ1]
			u	j J	*

где А – штрафной коэффициент; J* – множество индексов ограничений, за нарушение которых берется штраф.

Здесь также отметим, что если существует d D , при котором

max1 1 minu maxj J (ρзад − Bepξ2 [gj (d,u, ξ1, ξ2) ≤ 0])≤ 0,

ξ Ξ

то существует {d }≠ , при котором Ξ = Ξ и

∫ minu Mξ2 (•)P(ξ1)dξ1 = 0 .

Ξ\Ξ

При этом сформулированная задача (4.39), (4.40) переходит в двухэтапную задачу с жесткими ограничениями.

Задача 4. Имеются конструктивные и управляющие переменные. На этапе эксплуатации ХТП область неопределенных параметров та же, что и на этапе проектирования. Этот случай соответствует

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

Эта задача может быть сформулирована (в отличие от задачи 1) следующим образом):

C* = minMξ {C(d,u, ξ)}

d,u

при условии

max max gj (d,u, ξ) ≤ 0

ξ j

или

max gj (d,u, ξ) ≤ 0, j J.

ξ

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

Mξ {C(d,u, ξ)}≈ ∑νi C(d,u, ξ(i)),

l I1

где νi – весовые коэффициенты, ∑νi =1, I1 – множество аппроксимационных точек в области Ξ .

i I1

Совокупность точек ξ(i), i I1, будем обозначать через S1 , а множество критических точек на ν-м шаге – через S2(ν) = {ξ(k) :k I2(ν)}.

Алгоритм 3

Шаг 1. Положим ν =0. Выбираем совокупность аппроксимационных точек S1 и начальную совокупность критических точек S2(ν) .

Шаг 2. Решаем задачу

min ∑γiC(d,u, ξ(i)) ;

d,u i I1

gj (d,u, ξ(k)) ≤ 0, j =1, m; k I2(ν)

и определяем d(ν), u(ν) . Шаг 3. Решаем т-задач

max gj (d(ν),u(ν), ξk ), j =1, m

ξ Ξ

и определяем т точек ξ( j)*, j =1, m . Шаг 4. Образуем множество

R(ν) = {ξ( j)* : gj (d(ν),u(ν), ξ( j)*) > 0}.

Если это множество пустое, то решение задачи получено. В противном случае перейдем к шагу 5. Шаг 5. Определим

S2(ν+1) = S2(ν) R(ν).

Положим ν := ν +1 и переходим к шагу 2.

Характерной чертой алгоритма 3 является увеличение числа критических точек на каждом шаге, соответственно увеличивается число ограничений. Это является определенным недостатком, поскольку в некоторых случаях при большом числе критических точек число ограничений может стать слишком большим.

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

max gj (d(ν),u(ν), ξk ), j =1, m

ξ Ξ

находится в одной из вершин параллелепипеда Ξ [28]. В начальное множество критических точек S2(0) включается некоторое количество угловых точек куба Ξ , а на шаге 3 рассчитываются значения функций gj (d(ν),u(ν), ξk ), j =1, m во всех угловых точках куба Ξ , не принадлежащих множествам S2(ν) и S1 .

Среди этих точек выбираются m точек, в которых функции gj (d(ν),u(ν), ξk ),	j =		принимают наи-
		1, m
большие значения. Далее определим множество критических точек S2(ν+1) = S2(ν)	R(ν) и переходим к ша-
гу 2 алгоритма 2.

Задача 5. Имеются конструктивные и управляющие переменные. На этапе эксплуатации неопределенные параметры могут быть определены в каждый момент времени. Для обеспечения выполнения ограничений gj (d,u, ξ) ≤ 0, j J могут быть использованы конструктивные и управляющие переменные.

Для этого случая условие гибкости (работоспособности) можно записать в виде

ξ Ξ u, d j J [gj (d,u, ξ) ≤ 0]

или

χ(d) = max min max gj (d,u, ξ) ≤ 0.	(4.41)
ξ Ξ u j J

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

C* = Mξ minC(d,u, ξ)
d,u

при условии (4.41).

Используя квадратурную формулу, функцию Mξ {•} можно приближенно заменить выражением

Mξ {C*(d,u, ξ)}≈ ∑νiC*(d,u, ξi),

i I1

где νi – весовые коэффициенты; ξi – аппроксимационные точки; I1– множество индексов аппроксима-

ционных точек.

Операции суммирования и минимизации можно поменять местами и задача может быть представлена в виде

		C* = mini i	∑νiC(di,ui, ξi)
		d ,u	i I1
при выполнении условия гибкости и g	j	(di,ui, ξi) ≤ 0, i I , j J .
			1

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

Задача 6. Формулировка этой задачи та же, что и задачи 2, за исключением того, что условие гибкости (работоспособности) проекта ХТП записывается в жесткой форме

χ(d) = max min max gj (d,u, ξ) ≤ 0 .	(4.42)
ξ Ξ u j J

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

Использование возможности изменять параметры u (с помощью системы управления) на этапе функционирования процесса "облегчает" переменным d удовлетворять ограничениям, что, в свою очередь, позволит уменьшить коэффициенты запаса.

Двухэтапную задачу оптимального проектирования можно записать в виде

				gj (d,u, ξ) ≤ 0,
C* = min Mξ min C(d,u, ξ)					j J
d		u

при ограничениях (4.42).

Используя прием дискретизации, перепишем последнюю задачу в виде

		min(l)	∑νiC(d,u(i), ξ(i)),			(4.43)
		d,u	l I1
g		(d,u(i), ξ(i)) ≤ 0, j =			; i I ;	(4.44)
	j			1, m
			1
χ(d) = max min max gj (d,u, θ) ≤ 0 .						(4.45)
		ξ Ξ u U j J

Решение задачи (4.43) – (4.45) прямыми методами не представляется возможным, поскольку вычисление χ(d) в каждой точке может привести к очень большим объемам вычислений. В связи с этим здесь

будет рекомендована итерационная процедура, основанная на идеях метода "ветвей и границ" [44] и обеспечивающая приближение значений целевой функции (4.43). При этом не требуется непосредственно вычислять величину χ(d) .

В дальнейшем нам потребуются два соотношения:

min max f (x, y) ≥ max min f (x, y);				(4.46)
x	y	y	x
max max f (x, y) = max max f (x, y) ,				(4.47)
x	y	y	x

где x, y – векторы дискретных или непрерывных переменных.

Последнее соотношение является очевидным. Введем функцию

ϕ(d,u, ξ) = max gj (d,u, ξ).

j J

Тогда величина имеет вид

χ(d) = max min ϕ(d,u, ξ).

ξ Ξ u

В соответствии с соотношением (4.46) имеем

	χ(d) = maxmin ϕ(d,u, ξ) ≤ minmax ϕ(d,u, ξ) = χU (d),		(4.48)
	ξ Ξ u	u ξ Ξ
где χU = min max max gj (d,u, ξ) = min max max gj (d,u, ξ).
u ξ Ξ j J	u j J ξ Ξ
Введем обозначение
	ϕj(d,u) = maxgj (d,u, ξ),
	ξ
отсюда
	χU (d) = min max ϕj (d,u).
	u j
Известно, что задача вычисления χU (d) может быть сведена к следующей:
			(4.49)
	min α	ϕj (d,u) ≤ α.
	u,α
Из (4.48) следует, что если χU (d) ≤ 0, то χ(d) ≤ 0 . Поэтому условие
	χU (d) ≤ 0

является достаточным условием допустимости (работоспособности) проекта, определяемого вектором конструктивных параметров d. Для определения величины χU (d) необходимо решить задачу (4.49), то-

гда χU (d) = α* , где α * – оптимальное значение переменных α . Аналогично можно показать, что

χ(d) = max min max gj (d,u, ξ) ≥ χ(d),
ξ Ξ u U	j J
где
χL(d) = max max min gj (d,u, ξ) = max max min gj (d,u, ξ).
ξ Ξ j J u U	j J ξ Ξ u U

Отсюда следует, что если

χL(d) ≥ 0,

то и χ(d) ≥ 0 . Поэтому это условие является достаточным условием недопустимости (неработоспособности) проекта с вектором d. Для определения χ(d) необходимо решить т задач вида

max min gj (d,u, ξ), j				1, m	.
ξ Ξ u U
Каждая из этих задач эквивалентна следующей
max α			min gj (d,u, ξ) ≤ α.
ξ,α			u U
Таким образом, имеем
				χL (d) ≤ χ (d) ≤ χU (d) .		(4.50)

Следовательно, вычислив значения χL и χU , получим оценки снизу и сверху величины χ – критерия

гибкости (работоспособности) Гроссманна. Проанализируем физический смысл условия

χU (d) ≤ 0.

Будем искать такой вектор u , который обеспечивает допустимость вектора d при любых ξ :

u, ξ, j gj (d,u, ξ) ≤ 0.

Используя этот критерий, мы ищем единственный вектор u , который обеспечивает допустимость вектора d при любых значениях ξ . Напомним, что в критерии гибкости Гроссманна χ(d) каждому зна-

чению ξ соответствует свой вектор u , обеспечивающий допустимость вектора d.

Если разность χU − χL мала, то рассмотренный подход дает возможность оценить гибкость ХТП, в противном случае необходима какая-либо регулярная процедура, позволяющая изменить эту разность.

Рассмотрим одну из этих процедур [43, 44].
Разобьем область Ξ на N областей Ξi,	(i =		). Для каждой области определяем величину
		1, N
χUi	= min max max gj (d,u, ξ).
	u U j J			ξ Ξi
Для этого необходимо решить задачу
				min α	(4.51)
				u, α
	max gj (d,u, ξ) ≤ α.
	ξ Ξi
Определим теперь величину χU следующим образом:
χU = max min max max gj (d,u, ξ).
	i u U			j J ξ Ξi
Поскольку Ξi Ξ , то имеет место неравенство
откуда		χUi		≤ χU ,

χU = maxχUi ≤ χU . i

Далее можно показать, что

χ ≤ χU ≤ χU .

Следовательно, получена уточненная верхняя оценка критерия гибкости Гроссманна. Заметим, что чем плотнее покрытие области Ξ , тем ближе будет χU к χ. Однако такой путь может приводить к решению большого числа задач (4.51). В связи с этим рассмотрим другой путь вычисления χ(d). Для этого представим критерий χ(d) в виде

χ (d) = max ψ (d, ξ),

ξ Ξ

где ψ(d, ξ) = min max gj (d,u,ξ).

uj

Здесь вычисление χ(d) сводится к определению точки ξ* , в которой функция ψ(d, ξ) принимает

максимальное значение. Для определения этой точки воспользуемся процедурой метода "ветвей и границ" [48]. Цель этой процедуры будет состоять в том, чтобы разбивая область Ξ на все большее число подобластей Ξi , постараться локализовать точку ξ* .

Пусть на ν-м шаге область Ξ разбита на N областей Ξi(ν), i =1, Nν : Ξ = Ξ1(ν) Ξ(2ν) ... Ξ(Nνν) . Далее выбирается одна из областей Ξ(kνν) , которая в свою очередь разбивается на некоторое число областей. Для

простоты будем считать, что область Ξ(kνν) делится на две области: Ξ(Sν+1) и Ξ(qν+1) (Ξ(kνν) = Ξ(Sν+1) + Ξ(qν+1)) . В качестве области Ξ(kνν) берется та из областей Ξi(ν), (i =1, Nν ), в которой с наибольшей вероятностью нахо-

дится оптимальная точка ξ* .

Вычислим для каждой области Ξi(ν) величину χUi :

χUi ≥ max ψ (d, ξ).

ξ Ξi

Величина χUi является верхней оценкой для значения функции ψ(d, ξ) внутри области Ξi(ν) . Поэтому закономерно в качестве квазиоптимальной области выбрать область Ξ(kνν) , для которой величина χUi принимает наибольшее значение:

χUk	= maxχUi .
	i

Для проведения процедуры метода ветвей и границ на каждой итерации необходимо также знать нижнюю границу значения величины χ(d) = ψ(d, ξ*). Будем вычислять ее следующим образом [45]. Обо-

значим через ξ*i решение задач (4.51) и найдем

		ψ(d, ξ*i ) = min max gj (d,u, ξ*i ).
		u	j
Для этого необходимо решить задачу
						min α;		(4.52)
						u,α
		gj (d,u, ξ*i ) ≤ α, ( j =				).
					1, m
Вычислим ψ(d, ξ*i )	для всех областей Ξi(ν), (i =			).
			1, Nν
Введем величину
		χ(Lν) = max ψ (d, ξ*j ).
			j
Очевидно, что		χ(d) = max ψ(d, ξ) ≥ R(ν),
		ξ Ξ
что и определяет R (ν )	как нижнюю границу для максимального значения функции ψ(d, ξ). Пусть для
некоторой области выполняется		соотношение	R (ν ) ≥ χUl		(ν ) , тогда		в соответствии с	неравенством
χUl ≥ max ψ(d, ξ) имеем R(ν) ≥ ψ(d, ξ),		ξ Ξl(ν) . Следовательно, точка ξ*					заведомо не принадлежит области
ξ Ξi

Ξl(ν) и в дальнейшем не рассматривается. Процедура прекращается при выполнении соотношения

R(ν) −χUkν ≤ ε,

где ε – малая величина.

Если речь идет об оценке гибкости производства, а не о вычислении χ(d), то описанная процедура

может окончится раньше, чем выполнится последнее условие. Действительно, пусть на ν-й итерации выполнится условие

	max χi ≤ 0,
тогда	i
	χkν ≤ 0.
Далее, на каждом шаге необходимо найти два значения χ S		и χq , соответствующих областям Ξ(Sν+1)
и Ξ(ν+1)	, на которые разбивается квазиоптимальная область Ξ(ν)	. Для этого потребуется два раза решить
q	kν

задачу (4.51) и, кроме того, необходимо найти величины ψ(d, ξ*S ) и ψ(d, ξ*q) , дважды решив задачу (4.52) для i = S и i = q .

Вернемся теперь к решению задачи (4.43) – (4.45):

C*

= min M

minC(d,u, ξ)

g

j

(d,u, ξ) ≤ 0,

j J ;

(В)

B

ξ

d

u

χ(d) = max min max gj (d,u, ξ) ≤ 0.

ξ

u j

Используя полученные выше оценки χL (d),

χU (d) , можно получить оценки оптимального значения

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

C*

= min M

ξ

minC(d,u, ξ)

g

j

(d,u, ξ) ≤ 0,

j J ;

(Г)

Г

d

u

χU (d) ≤ 0.

CД* = min Mξ

{minC(d,u, ξ)

gj (d,u, ξ) ≤ 0,

j J };

(Д)

d

u

χL (d) ≤ 0.

Задачи (Г) и (Д) отличаются от задачи (В) только тем, что в них ограничение χ(d)≤0

заменено соот-

ветственно на ограничения χU (d) ≤ 0 и χL(d) ≤ 0. Поскольку имеет место неравенство

χL(d) ≤ χ(d) ≤ χU (d) ,

то можно записать

CД* ≤ CВ* ≤ CГ* ,

где CВ* ,CГ*,CД* – оптимальные значения целевой функции задач (В), (Г) и (Д), соответственно. Следует отметить, что решение задачи (Г) и (Д) проще, чем решение задачи (В). Если разность CГ* −CД* достаточ-

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

dk(B)* = 0,5(dk(Г)* + dk(Д)*)

при условии

χ(d(B)*) ≤ 0.

Введем еще одну вспомогательную задачу, разбив область Ξ на N областей Ξi (i =1, N ) и определяя

χUi	(d) = min max max gj (d,u, ξ).
	u U j J ξ Ξ
	CE* = min Mξ {minC(d,u, ξ)			gj (d,u, ξ) ≤ 0, j J };	(Е)
	d	u
	χUi (d) ≤ 0, ... , χUN (d) ≤ 0 .
Поскольку имеет место неравенство χ(d) ≤ χU ≤ χU (d) , то
	CB* ≤ CE* ≤ CГ* .
Пусть величина r(Ξi) характеризует размер подобласти Ξi . При выполнении условия
	r(Ξi) ≤ ε ,

где ε – достаточно малое число, можно получить достаточно хорошее приближение к решению задачи

(4.43) – (4.45).

Рассмотрим алгоритм решения задачи (4.43) – (4.45) с помощью задачи (Е), в которой разбиение на области Ξi будет проводиться более "экономичным" способом. Обозначим через Ξi(ν), i =1, N(ν) подобла-

сти, на которые разбивается область Ξ на k-й итерации.

Алгоритм 4 [44]

Шаг 1. Положим ν = 0. Выбрать начальное разбиение области Ξ на подобласти Ξi(ν), i =1, N(ν)) и начальное значение d (ν) вектора d .

Шаг 2. Решить задачу (Е). Пусть CE(ν) и d (ν) – оптимальные значения критерия и вектора d . Шаг 3. Найти множество S(ν) номеров активных ограничений:

χUi (d(ν)) = 0, i S(ν) .

Очевидны соотношения

χUi (d(ν) ) ≥ χUj (d(ν) ), i S(ν), j ≠ i .

Шаг 4. Если множество S(ν) – пустое, то решение задачи (4.43)-(4.45) получено. В противном случае перейти к шагу 5.

Шаг 5. Проверить условие

r(Ξi) ≤ δ, i S(ν) ,

где δ – заранее заданное малое число.

Если условие выполняется, то итерационную процедуру закончить, в противном случае перейти к шагу 6.

Шаг 6.		Разбить каждую область Ξi(ν) (i S(ν))		на две подобласти Ξi(ν+1)		и Ξi(ν+1)	и образовать новое
					1	2
разбиение, исключив из предыдущего разбиения подобласти					Ξi(ν) (i S(ν))	и добавив новые области
Ξi(ν+1)	, Ξi(ν+1)	(i S(ν) ).
1	2
Шаг 7. Положить ν:= ν+1 и перейти к шагу 2. Поскольку
		Ξi(ν+1)	Ξi(ν) , Ξi(ν+1) Ξi(ν) , χUi (ν)(d) ≥ χUi (ν+1)(d),		χUi (ν)(d) ≥ χUi (ν+1)(d) .
		1	2	1	2

Следовательно, CE(ν) ≥ CE(ν+1) .

Приведенный алгоритм позволяет определить локальный минимум задачи (4.43) – (4.45). Особенность этого алгоритма состоит в том, что на каждой итерации выполняется операция, кото-

рая приближает ограничение (Е) к ограничению (4.48) (шаги 5 и 6). Идея этой операции близка к идее метода "ветвей и границ" [48], поскольку на каждой итерации разбиению подвергаются те подобласти Ξi(ν) , для которых верна оценка величины χ(d) наибольшая. Фактически поиск можно прекратить при

выполнении условия

CE(ν) −CE(ν+1) ≤ ε,

где ε – достаточно малое число.

Задача 7. Формулировка этой задачи та же, что и задачи 3, за исключением того, что условие гибкости (работоспособности) проекта записывается в виде

χ(d) = min min max max gj (d,u, ξ) ≤ 0 .	(4.53)
ξ1 Ξ1 u U ξ2 Ξ2 j J

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

из области Ξ2 . Поэтому для фиксированного момента времени будем иметь следующую постановку оптимизационной задачи:

			C(d,u, ξ1, ξ2)
C(d, ξ1) = min M	ξ	2
u U

max gj (d,u, ξ1, ξ2) ≤ 0,	j J .
ξ2 Ξ2

В качестве критерия оптимального проектирования должно быть взято математическое ожидание

по ξ1 от величины C(d, ξ1) .

В результате приходим к задаче

C* = min M 1					(4.54)
			С(d, ξ1)
d	ξ

при ограничении (4.53).

Используя метод дискретизации критерия, получим дискретный аналог задачи (4.54), (4.53):

C* = mind,ui i∑I1 wilC(d,ui, ξ1i, ξ2l )

при условиях (4.53) и