авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 | 2 || 4 | 5 |   ...   | 8 |

«Д.С. ДВОРЕЦКИЙ, С.И. ДВОРЕЦКИЙ, Г.М. ОСТРОВСКИЙ НОВЫЕ ПОДХОДЫ К ПРОЕКТИРОВАНИЮ ХИМИКО-ТЕХНОЛОГИЧЕСКИХ ПРОЦЕССОВ, АППАРАТОВ И СИСТЕМ В УСЛОВИЯХ ...»

-- [ Страница 3 ] --

Поскольку математическое ожидание M { C (a, d, z, )} дает сред нее значение критерия интегрированного проектирования на стадии функционирования ХТС, то естественно использовать эту величину как целевую функцию задачи интегрированного проектирования в условиях неопределенности. Объединяя целевую функцию и условие гибкости ХТС – max g j (a, d, z, ) 0, j = 1, …, m, мы можем сформулировать ОЗИП с жесткими ограничениями в условиях неопределенности:

min M { С (a, d, z, )}, (2.1) a,d, z max g j (a, d, z, ) 0, j = 1,..., m. (2.2) Если функции распределения вероятностей для неизвестны, то можно использовать одну из следующих формулировок задачи:

I = min I (a, d, z ), a,d, z max g j (a, d, z, ) 0, j = 1,..., m, где I (a, d, z ) = wi С ( a, d, z, i ) или iI I (a, d, z ) = max С (a, d, z, ), (2.3) wi – весовые коэффициенты, i (i I1 ) аппроксимационные точки, I1 – множество индексов аппроксимационных точек.

ОЗИП (2.1), (2.2) имеет решение, если выполняется условие гиб кости проектируемой ХТС (a, d ) = min max max g j (a, d, z, ) 0. (2.4) a, d, z jJ Условие гибкости (2.4) гарантирует возможность удовлетворения всех ограничений (2.2) для всех значений из области.

При формулировании ОЗИП с мягкими (вероятностными) огра ничениями предположим, что мы имеем полную информацию относи тельно функции распределения вероятностей для. В этом случае ОЗИП имеет вид min I (a, d, z ), a,d, z { } P() d j, Pr g j (a, d, z, ) 0 = j = 1,..., m, (2.5) j j = { : g j (a, d, z, ) 0, }, где P() функция плотности вероятности.

В качестве целевой функции I(a, d, z) можно использовать (2.3), т.е. либо среднее значение целевой функции C(a, d, z, ) на стадии функционирования ХТС, либо наихудшее значение целевой функции интегрированного проектирования.

Главная трудность решения сформулированной ОЗИП состоит в необходимости вычисления многомерных интегралов M {С (a, d, z, )} P() d.

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

min, a, d, z, Pr { С (a, d, z, ) 0} 0, (2.6) { } Pr g j (a, d, z, ) 0 j, j = 1,..., m.

Здесь – новая переменная, ограничивающая целевую функцию С (a, d, z, ). Сравним формулировку ОЗИП (2.6) с постановкой ОЗИП (2.5). В задаче (2.6) мы ищем наименьшее значение переменной, для которой условие С (a, d, z, ) 0 удовлетворяется с заданной вероятностью 0. Таким образом, решив задачу ОЗИП (2.6), мы нахо дим конструкцию ХТС a, d и режим z, которые гарантируют, что в процессе функционирования ХТС целевая функция С (a, d, z, ) будет меньше, чем с вероятностью 0.

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

min, a, d, z, Pr { С (a, d, z, ) 0} 0, max g j (a, d, z, ) 0, j = 1,..., m.

Предположим, что мы решили ОЗИП и получили решение [ a, d, z ]. Для его реализации мы должны обеспечить выполнение условий z = z на стадии функционирования ХТС. Это означает, что мы лишены возможности подстраивать управляющие переменные z на стадии функционирования ХТС для выполнения регламентных требо ваний и проектных ограничений. Ясно, что постановка и решение од ностадийной задачи при проектировании приводят к не вполне эконо мичным конструкциям аппаратов ХТС, так как не допускается исполь зование настройки управляющих переменных z на стадии функциони рования ХТС.

Двухстадийная задача интегрированного проектирования.

Все формулировки двухстадийных задач интегрированного проекти рования (ДЗИП) будут учитывать возможность уточнения неопреде ленных параметров и настройки управляющих переменных z на ста дии функционирования ХТС. В этом состоит принципиальная разница между ДЗИП и ОЗИП [31]. В ОЗИП (см. задачи (2.1), (2.2) и (2.5)) пе ременные a, d, z равноправны в том смысле, что они не изменяются на стадии функционирования ХТС. В двухстадийной задаче возможны два случая: а) переменные a, d по-прежнему постоянны на этапе функ ционирования ХТС, в то время как переменные z могут изменяться;

б) переменная a и одна часть конструктивных переменных d k, k = 1, 2, …, k1 постоянны на стадии функционирования ХТС, в то вре мя как другая часть конструктивных переменных d k, k = k1, k1 + 1, …, K, и управляющие переменные z могут изменяться. В частности, это свойство позволяет настраивать конструктивные d k параметры наряду с управляющими z переменными для удовлетворения ограничений задачи.

При формулировке ДЗИП мы будем использовать следующее предположение – на стадии функционирования ХТС в каждый момент времени:

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

б) решается ДЗИП с использованием математической модели с уточненными неопределенными параметрами и найденный опти мальный режим реализуется на стадии функционирования ХТС.

Введем понятие области гибкости ХТС. Она состоит из точек об ласти неопределенности, для которых можно найти такие значения управляющих переменных z, при которых все ограничения задачи ин тегрированного проектирования g j (a, d, z, ) 0 будут выполняться.

Случай 1. На стадии функционирования ХТС определяются точные значения всех неопределенных параметров, при этом все ограничения являются жесткими.

Рассмотрим условие работоспособности (гибкости) ХТС: ХТС является гибкой, если для каждого можно найти такие значения (режимных) управляющих переменных z, что все ограничения задачи будут удовлетворены.

Задачу оптимизации в ДЗИП ХТС с использованием математиче ской модели y = (a, d, z, ) с уточненными неопределенными пара метрами на стадии функционирования ХТС будем называть внутрен ней задачей интегрированного проектирования. В данном случае она имеет вид С (a, d, ) = min C (a, d, z, ), (2.7) z g j (a, d, z, ) 0, j = 1,..., m.

Условие гибкости для задачи ДЗИП имеет вид 1 (a, d ) = max min max g j (a, d, z, ) 0.

jJ z Предположим, что функция плотности распределения вероятно сти P() известна. Поскольку в каждый момент времени на стадии функционирования ХТС значение критерия оптимизации будет равно С (a, d, ), то на стадии проектирования ХТС мы можем оценить бу дущую работу ХТС, подсчитав математическое ожидание M {•} вели чины С (a, d, ) :

{ } M C (a, d, ) = C (a, d, ) P() d.

Эта величина будет использоваться как целевая функция в задаче интегрированного проектирования в условиях неопределенности. Ре зультирующая двухстадийная задача есть задача стохастического про граммирования с рекурсией [32, 34]:

{ } min M C (a, d, ). (2.8) a,d Предположим, что внутренняя задача интегрированного проекти рования (2.7) имеет решение во всех точках и функция плотности распределения вероятности P() известна. Тогда имеем { } min min C (a, d, z, ) | g j (a, d, z, ) 0 P () d.

a,d z Поскольку интеграл есть бесконечная сумма и переменные z, со ответствующие различным, независимы друг от друга, то мы можем изменить порядок операторов интегрирования и минимизации:

min min C (a, d, z, ) P() d, (2.9) z () a,d g j (a, d, z (), ) 0, j = 1,..., m,.

Так как оптимальное значение z во внутренней задаче ДЗИП (2.7) зависит от, то z есть многомерная функция z(). Таким образом, в задаче (2.9) мы ищем оптимальные векторы a, d и оптимальную многомерную C (a, d, z (), ) P() d функцию z(), доставляющую функционалу минимальное значение. Объединяя оба оператора минимизации по a, d и z(), мы получим ДЗИП для случая 1 (ДЗИП1):

C (a, d, z(), ) P() d, min (2.10) a,d, z () g j (a, d, z (), ) 0, j = 1,..., m,. (2.11) Задача (2.10), (2.11) имеет бесконечное число ограничений и по исковых переменных (одна многомерная функция z() эквивалентна бесконечному числу обычных поисковых переменных). Решение зада чи (2.8) a*, d* гарантирует гибкость ХТС, так как внутренняя задача (2.7) ДЗИП1 решается во всех точках. При этом нельзя гарантировать, что задача (2.7) имеет решение для каждого a A, d D и.

Поэтому задача (2.8) должна быть дополнена условием гибкости 1 (a, d ) = max min max g j (a, d, z, ) 0, т.е.

jJ z {C (a, d, )}, C1 = min M a,d 1 (a, d ) 0.

В результате мы получаем другую постановку задачи ДЗИП1:

C (a, d, z (), ) P() d ), C1 = min a,d, z () g j (a, d, z (), ) 0, j = 1,..., m,, 1 (a, d ) = max min max g j (a, d, z, ) 0.

jJ z Заменим многомерный интеграл в целевой функции С1 некото рой конечной суммой с помощью соответствующей квадратурной формулы [35]:

C (a, d, z(), ) P() d = i i C (a, d, z, ), i i I где i (i I1 ) аппроксимационные (узловые) точки;

z i = z ( i ) век тор режимных (управляющих) переменных, соответствующий аппрок симационной точке i ;

i (i I1 ) весовые коэффициенты, удовлетво i = 1. Кроме того, заменим бесконечное ряющие условиям i 0 ;

iI число ограничений конечным числом ограничений только в аппрокси мационных точках i (i I1 ). Таким образом, получим дискретный вариант ДЗИП1:

i C ( a, d, z i, i ), C 1= min a, d, z i iI g (a, d, z i, i ) 0, i I1, 1 (a, d ) = max min max g j (a, d, z, ) 0.

jJ z В случае если функции распределения неопределенных парамет ров неизвестны, аппроксимационные точки желательно выбирать таким образом, чтобы они попадали в область наиболее вероятных значений, которые параметры могут принимать при функционировании ХТС, и достаточно плотно покрывали область. В некоторых случаях ДЗИП решается без ограничений 1 (a, d ) = max min max g j (a, d, z, ) 0 :

jJ z i C ( a, d, z i, i ), C1 = min a, d, z i iI g ( a, d, z i, i ) 0, i I 2.

Здесь множество I2 содержит аппроксимационные и некоторые дополнительные (критические) точки. Эта постановка задачи оправда на только в случае, когда аппроксимационные точки достаточно плот но покрывают область.

Это требует большого числа аппроксимационных точек даже для сравнительно малой размерности n вектора. Если число узловых точек по каждой компоненте вектора равно p, то число аппроксима n ционных точек будет равно p. В этом случае размерность задачи n будет равна na + nd + p n z. В случае когда число аппроксимационных точек невелико, использование ограничения 1 (a, d ) 0 совершенно необходимо, так как это гарантирует выполнение ограничений задачи не только в аппроксимационных точках i, i I1, но и во всех других точках области.

При интегрированном проектировании возможно оценить сред ние потери энергии, связанные с необходимостью выполнения регла ментных требований, проектных ограничений и неточности исходной математической модели ХТС – y = (a, d, z, ). Пусть a, d – реше ние ДЗИП1, а z (a, d, ) – решение внутренней задачи (2.7) при фиксированных a, d и параметре. Чтобы поддерживать значение управляющей переменной на уровне z (a, d, ), необходимо расхо довать энергию. Например, если z l – температура, то необходимо теп ло, чтобы поддерживать температуру zl ;

если z l – расход некоторого потока, то необходима энергия для насоса или компрессора, поддер живающего требуемое значение потока. Таким образом, величина z (a, d, ) непосредственно связана с потребляемой энергией.

Предположим, что потребляемая энергия пропорциональна величине z (a, d, ) с коэффициентом пропорциональности ki. Тогда среднее потребление энергии, связанное с реализацией оптимального значения i-й управляющей переменной, будет определяться величиной nz ki zi (a, d, ) P ( ) d.

I = i = В этом случае аналогично коэффициентам запаса конструктивных переменных можно ввести понятие энергетического коэффициента запаса E :

( ) E = I I N IN, где IN – потребляемая энергия при номинальных значениях неопреде ленных параметров.

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

В этом случае вектор неопределенных параметров состоит из двух подвекторов 1 и 2. Пусть при этом 1 и 22 ( = 1 U 2 ).

На стадии функционирования ХТС значения компонентов вектора могут быть измерены (определены) достаточно точно. В то же время компоненты вектора 2 не могут быть уточнены.

Условие гибкости для данного случая имеет вид [8] 2 (a, d ) = max min max max g j (a, d, z, 1, 2 ) 0, (2.12) 2 2 jJ 11 z причем 1 (a, d ) 2 (a, d ).

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

Эта ситуация соответствует ОЗИП с жесткими ограничениями. В этом случае формулировка внутренней задачи интегрированного проектирования имеет вид { } С (a, d, 1 ) = min E 2 С (a, d, z, 1, 2 ), z max g j (a, d, z, 1, 2 ) 0, j = 1,..., m, (2.13) 2 { } С ( a, d,, ) P ( ) d, E2 С (a, d, z, 1, 2 ) = 1 2 2 { } E2 С (a, d, z, 1, 2 ) где – математическое ожидание функции С (a, d, z, 1, 2 ) по переменной 2 при фиксированном 11. Условие 2 (a, d ) 0 гарантирует существование решения внутренней задачи интегрированного проектирования (2.13) для всех 11. Следователь но, оно должно быть использовано в постановке ДЗИП2 как ограниче ние. Так как величина { } E1 С (a, d, 1 ) = С (a, d, 1 ) P (1 ) d характеризует эффективность функционирования (будущую работу) ХТС, то она может служить в качестве целевой функции в ДЗИП2.

Сформулируем эту задачу в виде { }, C2 = min E1 С (a, d, 1 ) (2.14) a,d 2 (a, d ) 0. (2.15) Используя квадратурную для аппроксимации интеграла в целевой функции задачи и заменяя бесконечный набор ограничений конечным набором только в аппроксимационных точках, получим следующую задачу:

wi vk C (a, d, z i, 1,i, 2,k ), C2 = min (2.16) a, d, z i iI kK max g j (a, d, z i, 1,i, 2 ) 0, j = 1,..., m, i I1, (2.17) 2 2 (a, d ) 0, где 1,i, 2,k – аппроксимационные точки для внешнего и внутреннего интегралов в критерии оптимизации, z i = z (1,i ), wi (i I1 ) и vk (k K ) весовые коэффициенты, удовлетворяющие условию wi = 1;

v k = 1.

wi 0 ;

vk 0;

iI1 kK Если функция плотности распределения неизвестна, весовые ко эффициенты wi (i I1 ) и vk (k K ), а также множества аппроксима { } { } ционных S1 = 1,i : i I1 и критических S 2 = 2, k : k K точек должен назначать исследователь.

Следует заметить, что в ДЗИП1 вектор zi, соответствующий ап проксимационной точке i, используется только в этой точке (см. квад ратурную формулу (2.16)). С другой стороны, из (2.16), (2.17) видно, что в ДЗИП2 вектор zi, соответствующий точке 1,i, используется во всех аппроксимационных точках 1,i, 2,k из области { } (1,i ) = : 2,k 2 ;

1 = 1,i.

В работе [31] показано, что справедливо следующее неравенство:

C2 C1.

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

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

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

( ) min I1 (a, d ) = I1(1) (a, d ) + I1( 2) (a, d ), a,d P()d, 1 = : h(a, d, ) = min max g j (a, d, z, ) 0, jJ z C (a, d, ) = min C (a, d, z, ) | z d ) = 1 P() d, I11) (a, ( где g j (a, d, z, ) 0, j = 1, 2,..., m min C (a, d, z, ) P() d, I1 2 ) ( a, d ) = ( zZ \ \ 1 теоретико-множественная разность.

Недостаток сформулированной задачи состоит в том, что ограни чения g j (a, d, z, ) 0, j = 1, 2,..., m строго выполняются только для 1, а вне области 1 ограничения не должны выполняться. По этому можно получить решение ДЗИП, которому будет соответство вать большие нарушения ограничений. Более того, ДЗИП может вооб ще не иметь решения. Поэтому было бы лучше, если нарушение этих ограничений было распределено равномерно по всей области.

В связи с этим дадим другую формулировку задачи оптимизации с мягкими ограничениями.

Введем штрафную функцию m p( g j (a, d, z, )), С (a, d,, ) = С (a, d, z, ) + j = 0, если g j 0, где p( g j ) = 2 штрафной коэффициент.

g j, если g j 0, Далее мы будем использовать следующую внутреннюю задачу оптимизации для всех значений :

С (a, d,, ) = min С (a, d,, ).

z Ясно, что оптимальное значение целевой функции этой задачи за висит от a, d,,, и пусть z(a, d,, ) – решение этой задачи. Тогда пе ременные a, d, должны быть выбраны таким образом, чтобы ограни чения удовлетворялись с заданной вероятностью и среднее значение величины C (a, d,, ) принимало минимальное значение. ДЗИП в этом случае будет иметь вид { } I = min E С (a, d,, ), a,d, P() d, 2 = { : g j (a, d, z (a, d,, ), ) 0}.

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

Случай 4. На стадии функционирования ХТС можно опреде лить точные значения неопределенных параметров, при этом имеются смешанные ограничения: ограничения с номерами j J1 = {,..., m1} являются жесткими, а ограничения с номерами j J 2 = {m1 + 1,..., m} – мягкими и должны быть удовлетворены с заданной вероятностью.

В этом случае ДЗИП будет иметь вид min I (a, d ) = I1 (a, d ) + I 2 (a, d ), a,d Pr[ 1 ], 1 (a, d ;

J1 ) = max min max g j (a, d, z, ) 0, jJ z где С I1 ( a, d ) = (a, d, ) P()d, С I 2 ( a, d ) = (a, d, ;

J 1 ) P() d;

\ С (a, d, ) = min С (a, d, z, ), С (a, d, ) решение задачи z g j (a, d, z, ) 0, j = 1,..., m;

и С (a, d, ;

J1 ) = min С (a, d, z, ), С (a, d, ;

J 1 ) решение задачи z g j (a, d, z, ) 0, j J 1.

Заметим, что задача с критерием I1(•) решается только для 1, а задача с критерием I2(•) – в остающейся части области. Ограниче ние Pr [ 1 ] гарантирует удовлетворение всех ограничений с вероятностью. Кроме того, дополнительное ограничение 1 (a, d ;

J1 ) = max min max g j (a, d, z, ) 0 гарантирует жесткое удов jJ z летворение ограничений с номерами j = 1, …, m1 для всех точек области.

2.2. ВЫЧИСЛЕНИЕ И АНАЛИЗ СВОЙСТВ ФУНКЦИИ ГИБКОСТИ ХТС Функция гибкости ХТС позволяет оценить ее работоспособность и является важнейшей характеристикой ХТС. Вычисление ее значений связано с решением сложной оптимизационной задачи. Для оценки гибкости ХТС достаточно знать только знак функции гибкости, поэто му далее будем различать следующие две задачи: 1) определение знака функции гибкости;

2) вычисление значений функции гибкости.

Для вычисления значений функций гибкости 1 (a, d ) и 2 (a, d ) будем применять методы верхней и нижней оценок функций гибкости с использованием схемы метода ветвей и границ, описанной в книге [31]. Метод ветвей и границ является итерационной процедурой [36, 37], которая, разбивая область на подобласти l, позволяет определить положение глобального минимума задачи:

I = min C (), g j () 0, j = 1,..., m, (2.18) { }, = i : iL i U i где I глобальный минимум функции C () в допустимой области { } D = :, g j () 0, j = 1,..., m задачи (2.18), функции g j () – непрерывно дифференцируемы.

Далее запишем задачу (2.18) для некоторой подобласти { } l = i : iL,l i U,l, i = 1,..., n, i I l = min C (), l (2.19) g j () 0, j = 1,..., m, где I l глобальный минимум функции C () в допустимой области { } Dl = : l, g j () 0, j = 1,..., m.

Будем называть величину µ l нижней оценкой (границей) опти мального значения I l целевой функции задачи (2.19), если выполня ется неравенство µ l I l (2.20) или µ l C () для всех l.

Желательно, чтобы нижняя оценка обладала следующими свойст вами.

Свойство 1. Нижняя оценка µl должна быть близка к I l.

l q, µl µ q, Свойство 2. Если то и если l = l1 U l2 U... U l N, то li, µ li µ l, для всех l li. Другими словами, здесь требуется, чтобы разбиение приводило к улучшению нижней оценки.

Свойство 3. Должно выполняться следующее соотношение:

lim µ l = C ( l ), (2.21) r ( l ) lim l = l, если r ( l ) где l точка, принадлежащая l, в которую стягивается l при r ( l ) 0, а r ( l ) величина, характеризующая размер подобласти l. Например, r ( l ) может иметь вид r ( l ) = max (U,l iL,l ).

i Свойство 4. Трудоемкость вычисления нижней оценки должна быть существенно меньше трудоемкости решения задачи глобальной оптимизации (2.18).

Аналогично можно ввести понятие верхней оценки l величины I l :

l I l. (2.22) Следует заметить, что оценки µl и l дают глобальную инфор мацию относительно функции I l в области l.

Предположим, что на -й итерации область разбита на N по добластей l (l = 1,..., N ). Обозначим множество подобластей l через ( ). Какие подобласти должны выбираться для разбиения? Как уменьшить число подобластей, которые должны рассматриваться?

Большинство процедур метода ветвей и границ основывается на следующих правилах.

1. На -й итерации область p исключается из рассмотрения, ес l, что выполняется соотноше ли найдется такая подобласть ние µ p l. В этом случае не имеет смысла сохранять область p, поскольку глобальный минимум заведомо не принадлежит этой облас ти.

2. На -й итерации подобласть l исключается из рассмотрения, если выполняется условие j l g j () 0. В этом случае вся подобласть l принадлежит недопустимой области задачи (2.18).

3. На -й итерации разбивается область lk, в которой нижняя граница принимает наименьшее значение, т.е. имеет место соотноше ние [36 – 38]:

µ l = min µ l, (2.23) lL где L = (1,..., N ) множество номеров подобластей l на -й ите рации.

Правило 3 основывается на следующем предположении: если µ p µ l, то I I l. Это правило является эвристическим. Ясно, что p если нижняя оценка близка к минимальному значению целевой функ ции задачи (2.19) во всех подобластях l, то последнее предположе ние будет выполняться в большинстве случаев. Однако, если нижняя оценка является грубой, то предположение может нарушиться. В этом случае использование этого правила приводит к дроблению подобла сти l, не содержащей глобальный минимум. Поэтому все подобласти за исключением подобласти l и подобластей, для которых условие µ p l выполняется, должны сохраняться как потенциальные канди даты, которые, вообще говоря, могут содержать глобальный минимум.

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

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

Вторая группа состоит из подобластей, для которых выполняется условие µ p l. Эти области могут быть удалены из дальнейшего рассмотрения.

Третья группа содержит оставшиеся подобласти.

Наименьшая из верхних оценок на -й итерации находится из ре шения задачи = min l, (2.24) lL где L = (1,..., N ).

Основанный на приведенных выше правилах алгоритм метода ветвей и границ имеет следующий вид [31].

Алгоритм 2. Шаг 1. Задаем i (i = 1,..., n),, начальное разбиение (1) на по добласти i (i = 1,..., N1 ). В частности, N1 может быть равно 1.

В этом случае (1) = (1) = (). Полагаем = 1.

Шаг 2. Рассчитываем нижние µl и верхние l (l = 1,..., N1 ) оцен ки всех подобластей i (i = 1,..., N1 ).

Шаг 3. Определим подобласть с наименьшим значением нижней оценки. Номер этой подобласти равен l. Таким образом, µ l µ l для всех l l.

Шаг 4. Прекращаем поиск, если выполняется одно из следующих условий:

1) U,l iL,l i, i = 1,..., n ;

i 2) l µ l. (2.25) Шаг 5. Разбиваем подобласть l на две подобласти p и q ( l = p U q ), ввести переменную s и точку cs ветвления:

{ } p = : l, s c s, = { : l, s cs }.

q Можно использовать различные способы выбора переменной ветвления и точки ветвления. В частности, в качестве переменных ветвления можно поочередно использовать все переменные или вы бирать переменную, соответствующую наибольшему ребру подобла сти l.

Шаг 6. Определим нижние µ p, µ q и верхние оценки p, q для подобластей p и q.

Шаг 7. Рассчитываем = min ( 1,, ).

p q Шаг 8. Проверяем выполнение условия µ j, j = p, q. По добласть, для которой это условие выполняется, удаляем из рассмот рения.

Шаг 9. Образуем множество +1. Для этого удаляем l из множества подобластей p и q, если они не были удалены на шаге 8.

Шаг 10. Полагаем := + 1 и переходим к шагу 3.

Сделаем некоторые пояснения.

На шаге 4 условие 2) гарантирует, что с точностью до глобаль ный минимум С () находится в подобласти l. Действительно, пусть выполняется условие µ l µ l для всех l l. Из этого усло вия и условия (2.25) вытекает следующее условие: µ l I l для всех l l, а также условие I l I l для всех l l. Таким образом, глобальный минимум C () находится в подобласти l, причем с точностью до он равен µ l.

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

На -й итерации возможны следующие два случая.

1. µ r µ l для всех l l, где µ r = min (µ p, µ q ). (2.26) На ( + 1) -й итерации в соответствии с алгоритмом метода ветвей и границ будет разбиваться подобласть l +1 (l +1 = r ) (полученная разбиением подобласти lk на данной итерации). Это желательный результат, и в этом случае будем говорить, что итерация является ус пешной. Условие (2.26) будет выполняться, если в каждой подобласти нижняя оценка совпадает с глобальным минимумом. Поэтому итера ция всегда будет успешной, если в каждой подобласти мы можем най ти такую нижнюю оценку. В этом случае алгоритм метода ветвей и границ будет всегда разбивать подобласть, содержащую глобальный минимум.

2. µ r min µ l. В этом случае на ( + 1) -й итерации в соответ l : l l ствии с алгоритмом будет разбиваться подобласть, которая была обра зована на одной из предыдущих итераций. Это так называемый «об ратный ход» алгоритма. Обратный ход будет встречаться тем чаще, чем более грубо с помощью нижней оценки оценивается глобальный минимум задачи (2.19). В худшем случае алгоритм метода ветвей и границ может превратиться в процедуру простого перебора.

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

Пусть в ходе выполнения алгоритма метода ветвей и границ было сделано итераций. Процедуру метода до -й может быть представ лена в виде дерева-графа (рис. 2.1.), у которого каждая вершина Ai j соответствует некоторой подобласти ij области.

Здесь верхний индекс j есть номер итерации (при j 1 ), на ко торой была получена эта вершина (подобласть), а нижний индекс есть номер подобласти, полученной на j-й итерации (в отличие от преды дущего текста на рисунке подобласти обозначаются через ij ). Об ласти 0 ( 1 = ) соответствует вершина A10. Вершину, не имеющую i потомков, будем называть висячей (на рис. 2.1 это вершины ( A12, A1, A2, A14, A2 ). Каждая вершина Ai j является либо висячей, 3 3 либо имеет два потомка A1, A2, соответствующих подобластям, кото рые были получены из области ij, соответствующей самой вершине Ai j. Если вершины A1k, A2 являются потомками вершины Ai j и полу k А А А А А А14 А А1 А Рис. 2.1. Представление метода ветвей и границ в виде дерева-графа чены на следующей итерации, то = j + 1. Нижнюю оценку, соответ ствующую области ij, мы будем обозначать через µ ij.

На рисунке 2.1 представлено дерево-граф, соответствующее че тырем шагам метода ветвей и границ. Пусть на первом шаге первой итерации имеется только одна подобласть 1, совпадающая с (рис. 2.2, а), и пусть на первой итерации область 1 разбита на подоб ласти 1, 1 (рис. 2.2, б). Вершины A1, A2 соответствуют подобла 1 1 стям 1 и 1.

1 Предположим, мы вычислили нижние оценки µ1, µ1 для этих 1 подобластей. Пусть µ1 µ1. Тогда согласно стратегии метода ветвей и 2 границ (шаги 3, 4 и 5 алгоритма 2.1) на второй итерации мы должны разбить область 1 на две подобласти 1 и 2. Таким образом, после 2 второй итерации имеется три подобласти 1, 2, 3 (рис. 2.2, в), кото 2 рым соответствуют вершины A12, A2, A1 (причем 3 = 1 ).

2 1 Пусть выполняются соотношения µ 2 µ1, µ 2 µ 3. Тогда на 2 2 третьей итерации мы должны разбить подобласть 2 на подобласти 3 и 3.

2 А10 А1 1 А А1 А1 А 2 3 А13 А14 5 3 А 4 А 2 1 1 1 1 1 3 0 2 1 2 3 А2 3 А2 А2 6 4 А 2 2 4 4 в) а) б) г) д) Рис. 2.2. Разбиение области на подобласти на итерациях 1, 2, 3, Таким образом, после третьей итерации область будет разбита на 4 подобласти: 1, 3, 3, 3 (рис. 2.2, г), которым соответствуют 2 3 вершины A1, A12, A1, A2. Область 3 соответствует вершине A2.

1 3 3 Пусть µ1 µ 3, i = 2, 3, 4, тогда на 4-й итерации должна быть разбита i область 1 на подобласти 4, 5, т.е. вершина A1 будет иметь два 3 4 потомка A14, A2. Таким образом, после 4-й итерации мы будем иметь подобластей i4, i = 2, …, 6 (рис. 2.2, д), которым отвечают вершины A12, A1, A2, A14, A2, соответственно.

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

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

~ Функцию L (C ;

) будем называть нижней оценочной функцией C в области, если во всех точках области выполняется неравен ~ ~ ство L (C ;

) C (). Если нижняя оценочная функция L (C ;

) являет ся выпуклой функцией, будем называть ее выпуклой нижней оценоч ной функцией и обозначать через L(C;

). Требуется, чтобы выпуклая нижняя оценочная функция удовлетворяла следующим условиям:

1) L(C ;

) должна быть дважды дифференцируемой функцией.

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

2) должно удовлетворяться следующее условие:

lim max | L(C ;

) – С| = 0.

r ( ) Аналогично можно ввести понятие вогнутой верхней оценочной функции U (C;

), удовлетворяющей условиям:

1) U (C;

) C () для всех ;

2) U (C ;

) вогнутая функция в области.

Выпуклые нижние оценочные функции могут быть использованы для вычисления нижних оценок. Действительно, рассмотрим задачу µ l = min L[C ();

l ], (2.27) l [ ] L g j ();

l 0, j = 1,..., m.

l решение этой задачи. Поскольку функции Пусть [ ] L[C ();

l ], L g j ();

l ( j = 1,..., m) выпуклы, то локальный минимум задачи (2.27) совпадает с глобальным минимумом этой задачи. Поэто му решение может быть получено с использованием обычных методов нелинейного программирования. Поскольку выполняются условия [ ] L[C ();

l ] C (), L g j ();

l g j (), j = 1,..., m для всех l, то решение задачи (2.27) дает нижнюю оценку оптимального значения критерия задачи (2.19), т.е. µ l l. При этом нижняя оценка не стано вится хуже при дроблении области (обычно она улучшается) [31].

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

Пусть l k -я подобласть разбивается на p -ю и q -ю подобласти и µ p µ q. В соответствии с шагом 3 для разбиения мы должны выбрать подобласть q. Это разбиение будет успешным, если I I q. Поэто p му желательно выбрать в качестве переменной разбиения переменную, для которой (µ p µ q ) является наибольшей величиной. Однако до разбиения нам неизвестны величины µ p и µ q. Поэтому необходимо выбрать косвенную оценку величины (µ p µ q ). Такой оценкой может быть чувствительность целевой функции С () по отношению к пере менным i в точке l решения задачи (2.27).

Функция гибкости 1 (a, d ) может быть представлена в виде 1 (a, d ) = max h(a, d, ), h(a, d, ) = min max g j (a, d, z, ).

где (2.28) zZ jJ Таким образом, вычисление функции гибкости сводится к макси мизации функции h(a, d, ) по переменным. Поскольку метод вет вей и границ изложен применительно к задаче минимизации, мы будем рассматривать задачу поиска минимума функции I = h ( a, d, ).

Следуя книге [31], рассмотрим методы вычисления верхней и нижней оценок функции гибкости 1 (a, d ). Изменяя порядок выпол нения операторов max и min, получим некоторую функцию 1 (a, d ) = min max max g j (a, d, z, ).

U Согласно равенству zZ jJ max max f ( x, y ) = max max f ( x, y ) = max f ( x, y ) операторы max и x y y x x, y jJ max перестановочны, поэтому имеем:

1 (a, d ) = min max max g j (a, d, z, ).

U (2.29) zZ jJ Свойство 1. Функция 1 (a, d ) определяет некоторую верхнюю U границу функции 1 (a, d ) в области :

1 (a, d ) 1 (a, d ).

U (2.30) Это свойство непосредственно следует из неравенства min max f ( x, y ) max min f ( x, y ). Из соотношений (2.28) и (2.30) по x x y y лучаем h(a, d, ) (a, d ),. Так как последнее неравенство U справедливо для любой области, то оно справедливо для любой по добласти i. Отсюда получаем h(a, d, ) 1,i (a, d ), 1, U (2.31) где 1,i (a, d ) = min max max g j (a, d, z, ), U zZ jJ i 1,i (a, d ) = max h(a, d, ) 1,i (a, d ).

U и (2.32) i Свойство 2. Пусть область разбита на подобластей Nk i( k ) (i = 1,..., N k ). Справедливы следующие неравенства [31]:

1,( k ) (a, d ) 1 (a, d ), U U (2.33) 1,( k ) (a, d ) 1 (a, d ), U (2.34) где 1,( k ) (a, d ) = max 1,i (a, d ) = max min max max g j (a, d, z i, ), (2.35) U U iL iL jJ i zi L = (1,..., N k ), k номер итерации в методе ветвей и границ.

Свойство 3. Величина 1,( k ) (a, d ) может быть получена реше U нием задачи 1,i (a, d ) = min u, U z,u max g j (a, d, z, ) u, j = 1,..., m. (2.36) i [z, u ] решение g j ( a, d, z ) = i i Пусть этой задачи и = max g j (a, d, z, ). Тогда выражение для 1,i (a, d ) U может быть пере i писано в виде 1,i (a, d ) = min max g j (a, d, z ), что эквивалентно сле U jJ z дующей задаче:

1,i (a, d ) = min u, U z,u g j (a, d, z ) u, j = 1,..., m.

Свойство 4. Имеет место следующее соотношение:

lim 1,i (a, d ) = h(a, d, i ), U (2.37) r ( i ) где r ( i(k ) ) некоторая мера размера подобласти i(k ) ;

i такая точ i(k ), i(k ) ка подобласти в которую стягивается при r ( i( k ) ) 0 ( i( k ) ). i Равенство (2.37) следует непосредственно из вида функции 1,i (a, d ). Таким образом, если подобласть i достаточно мала U (т.е. r ( i ), где достаточно малая величина), то 1,i h(a, d, i ), i i, U (2.38) h(a, d, i ) = min u, z,u g j (a, d, z, ) u, j = 1,..., m.

i (2.39) Поскольку функция 1,i (a, d ) есть решение задачи U 1,i (a, d ) = min u, U z max g j (a, d, z, ) = u, j J A (a, d ;

i ), i то функция h(a, d, i ) определяется только множеством J A (a, d ;

i ) активных ограничений задачи (2.39).

Свойство 5. Пусть область разбита на N k подобластей i( k ) (i = 1,..., N k ). Обозначим через (k ) множество подобластей i( k ) (i = 1,..., N k ) и через (k ) некоторое подмножество множества (k ), которое содержит подобласть с наибольшей верхней границей 1,i (a, d ). Пусть на k-й итерации разбиваются все подобласти под U множества (k ). В результате получаем подмножество (k ). Таким образом, на (k + 1)-й итерации мы получим новое подмножество по добластей:

{ } ( U ( ( k ) / ( k ) ) ).

( k +1) = ( k +1), i = 1,..., N k +1 ( k +1) = ( k ) i Тогда выполняется следующее неравенство:

1,( k ) 1,( k +1).

U U (2.40) Действительно, новое множество подобластей i( k +1) обладает следующим свойством:

i( k +1) ( k ) (jk ) ( k ) : i( k +1) (jk ).

Поскольку i( k +1) i( k ), то из вида функции 1,i (a, d ) следует U неравенство 1,i (a, d ) 1, j. Поскольку подмножество (k ) содержит U U подобласть с наибольшей верхней границей 1,i, то получим U 1,( k ) = max 1, j (a, d ) max 1,i (a, d ) = 1,( k +1).

U U U U j i На практике соотношение (2.40) за исключением особых случаев является строгим неравенством, поэтому дробление таких подмно жеств улучшает верхнюю границу 1,( k ).

U Рассмотрим теперь метод вычисления нижней границы функции 1,i (a, d ) в подобласти i.Запишем следующее неравенство 1,i (a, d ) = max h(a, d, ) max h(a, d, l ), i lS i где S i произвольное множество точек из подобласти i ( S i i ).

1,i (a, d ) = max h(a, d, l ) есть нижняя L Таким образом, величина lS i граница функции 1,i (a, d ) для любого множества S i. Правая часть этого равенства – решение задачи min, max h(a, d, l ), lS i которую можно путем последовательных преобразований свести к окончательному виду задачи вычисления нижней границы функции 1,i (a, d ) :

1,i (a, d ) = min, L z l, g j (a, d, z l, l ), j = 1,..., m, l S i, (2.41) где z l новые поисковые переменные.

Множество S i желательно выбирать таким образом, чтобы ниж няя граница 1,i (a, d ) была по возможности ближе к величине L h(a, d, i ), где точка i точка максимума функции h(a, d, ) по пе ременной в подобласти i.

Отметим следующие факты:

1) в соответствии с методом ветвей и границ дробятся подобла сти, которые могут содержать максимум функции h(a, d, ), и поэто му теоретически размер подобласти, содержащей максимум функции, должен стремиться к нулю;

2) величина h(a, d, i ) определяется только множеством актив i ных точек S AP ;

J A ( a, d ;

i ) J A (a, d ;

i ), если 3) множество стремится к r ( i ) 0 и lim i = i.

r ( i ) Поэтому в качестве множества S i разумно использовать множе i ство S AP, которое состоит из точек, соответствующих активным огра ничениям задачи (2.36).

Поскольку для каждой подобласти верно неравенство 1,i (a, d ) 1,i (a, d ), то мы получаем L max 1,i (a, d ) max 1,i (a, d ) 1 (a, d ).

L i i Таким образом, величина 1,( k ) (a, d ) = max 1,i L L (2.42) i является нижней границей функции 1 (a, d ).

Функция гибкости 2 (a, d ) может быть представлена в виде 2 (a, d ) = max h ( 2) (a, d, 1 ), (2.43) где h ( 2) (a, d, 1 ) = min max max g j (a, d, 1, 2 ). Вычисление функции jJ 2 z гибкости 2 (a, d ) сводится к максимизации функции h ( 2) (a, d, 1 ) в области 1 по переменным 1. Сведем задачу (2.43) к виду h ( 2) (a, d, 1 ) = min, (2.44) z, max g j (a, d, z, 1, 2 ), j = 1,..., m.

2 Подобно функции 2 (a, d ) можно ввести функцию гибкости 2,i (a, d ) для подобласти 1 области i 2,i (a, d ) = max min max max g j (a, d, z, 1, 2 ). (2.45) jJ 11 2 z i Аналогично формуле (2.43) эта функция может быть представле на в виде 2,i (a, d ) = max h ( 2) (a, d, 1 ). (2.46) Введем вспомогательные функции U (a, d ) и U,i (a, d ), соот 2 ветствующие всей области 1 и любой ее подобласти 1 :

i U (a, d ) = min max max max g j (a, d, z, ), 11 2 2 jJ z U,i (a, d ) = min max max max g j (a, d, z, ).

11 2 2 jJ z i Эти функции получены изменением порядка выполнения двух операторов в функциях 2 (a, d ) и 2,i (a, d ). Свойства функций U (a, d ) и U,i (a, d ) аналогичны свойствам функций 1 (a, d ) и U 2 1,i (a, d ) [31].

U Рассмотрим теперь проблему вычисления нижней границы функ ции гибкости 2 (a, d ) в некоторой подобласти 1 1. Используя те i же самые рассуждения, которые мы применяли для определения ниж ней границы функции гибкости 1 (a, d ), мы получим следующую за дачу для определения нижней границы функции гибкости 2 (a, d ) в подобласти 1 :

i 2,i (a, d ) = min u, L z l,u max g j (a, d, z l, 1,l, 2 ) u, j = 1,..., m, 1,l S i, 2 где S i произвольное множество точек 1,l в области 1. Разумно брать точки j,i ( j J A (a, d ;

1 )) в качестве множества S i в сформу i лированной выше задаче. Имея алгоритмы для вычисления верхней и нижней границы функции гибкости 2 (a, d ), мы можем вычислять значения функции гибкости, применяя метод ветвей и границ для мак симизации функции h( 2) (a, d, 1 ) в области 1 по переменным 1.

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

Будем предполагать, что выполняется следующее условие: функ ции g j (a, d, z l, lj ) 0 для всех j, l являются строго квазивыпук лыми по переменным z. Для вычисления функции гибкости будем использовать метод ветвей и границ со следующими изменениями.

Вместо задачи (2.36) для вычисления верхней границы функции h(a, d, ) мы будем решать следующую задачу [6]:

1,i = min u, U zZ,u max U ( g i ;

i ) u, j = 1,..., m, (2.47) i где U ( g i ;

i ) является строго квазивыпуклой функцией по перемен ным z и верхней оценочной вогнутой функцией для g j (a, d, z, ) по переменным, удовлетворяющей условиям U ( gi ;

i ) g j, i, [ ] lim max U ( g j ;

i ) g j (a, d, z, ) = 0, z Z, a A, d D. (2.48) r ( i )0 i Для определения нижней границы будем решать задачу (2.41), локальный минимум которой совпадает с глобальным минимумом.

Сравнивая задачу (2.36) и (2.47) и учитывая условие (2.48), получим следующее неравенство: 1,i 1,i h(a, d, ), i. Из условия U U (2.48) и задачи (2.47) следует, что имеет место следующее равенство:

lim 1,i = h(a, d, i ), U (2.49) r ( i ) где i стягивается в точку i. Это означает, что, когда подобласть i уменьшается, величина 1,i становится все более точной верхней U h ( a, d, ).

оценкой значения Из следует равенство (2.36) lim 1,i = h(a, d, i ), где i i при r ( i ) 0. Следовательно, U r ( i ) получаем lim 1,i = lim 1,i = h(a, d, i ` ).

U U r ( i )0 r ( i ) Обозначим через i(k ) подобласть, имеющую наибольшую верх нюю границу на k-й итерации. Так как на каждой итерации в соответ ствии с методом ветвей и границ дробится подобласть i(k ), то при r ( i( k ) ) 0 значение 1,i стремится к 1,i. Следовательно, в пределе в U U задаче (2.47) мы будем получать то же решение, что и в задаче (2.36).

Таким образом, метод ветвей и границ, использующий в качестве нижней и верхней границы решения задач (2.41) и (2.47) позволяет получить глобальное решение задачи вычисления значения функции гибкости 1 (a, d ).

Описанная модификация метода ветвей и границ для определения глобального решения задачи вычисления значения функции гибкости 1 (a, d ) сравнительно легко обобщается для вычисления функции гибкости 2 (a, d ). Однако для этого необходимо, чтобы функции g j (a, d, z, ) ( j = 1,..., m) были строго квазивогнутыми по параметрам 2.

2.3. АЛГОРИТМЫ РЕШЕНИЯ ОДНО- И ДВУХСТАДИЙНЫХ ЗАДАЧ ИНТЕГРИРОВАННОГО ПРОЕКТИРОВАНИЯ ХТС С «ЖЕСТКИМИ» ОГРАНИЧЕНИЯМИ Как и ранее, здесь будем использовать следующие обозначения:

{ }, = L U, – вектор неопределенных парамет ров, принадлежащих области ;

причем = {1, 2}, где 1 1 – под вектор компонентов, которые могут быть с достаточной точностью определены на стадии эксплуатации химико (измерены) технологической системы;

2 – подвектор компонентов, которые не удается идентифицировать даже на стадии эксплуатации ХТС;

{ } = i i = 1, – множество производимых продуктов (ассортимент);

d D – вектор проектных (конструктивных) параметров (множество D определяется типом a A аппаратурного оформления химического производства);

z Z – вектор режимных (управляющих) переменных ХТС;

C(•) – критерий оптимального проектирования ХТС.

Перейдем к рассмотрению задачи интегрированного проектиро вания ХТС, в которой конструктивные переменные a A, d D и ре жимные (управляющие) переменные z Z должны быть выбраны та ким образом, чтобы минимизировать приведенные затраты C(•), вклю чающие стоимость реализуемого проекта (капитальные затраты) и за траты на эксплуатацию ХТС.

Эксплуатационные затраты включают в себя: 1) стоимость сырья и материалов;

2) стоимость потребляемой технологическим оборудо ванием электро- и тепловой энергии;

3) заработную плату обслужи вающего персонала;

4) затраты на социальные нужды;

5) затраты на содержание и эксплуатацию технологического оборудования.

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

1) убедиться в гибкости проекта при найденных векторах a*, d* в задаче интегрированного проектирования, т.е. показать, что 1 ( a, d ) 0 ;

2) максимизировать меру гибкости многоассортиментного про изводства и в то же время минимизировать стоимость проекта.

Формулировку условия гибкости (задающего ограничения задачи) определяют следующие факторы.

1. Характер информации, содержащей неопределенность. Неоп ределенность может быть параметрической или модельной. В первом случае известна форма математической модели, но неизвестны точные значения некоторых ее параметров. Во втором случае предполагают, что нет точного знания о «модели» (триаде «модель–алгоритм– программа») ХТС. Имеется ряд альтернативных «моделей» (триад «модель–алгоритм–программа»), одна из которых соответствует дей ствительности.

2. Существование и величина неопределенности информации на второй стадии (на первой стадии неопределенность присутствует практически всегда).

3. Способ обеспечения гибкости ХТС:

• имеются конструктивные и управляющие переменные;

• имеются только конструктивные переменные;

• имеются только управляющие переменные.

4. Тип ограничений: ограничения могут быть «жесткими» (без условными) и «мягкими» (вероятностными). Жесткие ограничения не должны нарушаться ни при каких условиях. Мягкие ограничения должны выполняться с заданной вероятностью. В нашей работе мы будем рассматривать следующие случаи:

• все ограничения являются «жесткими»;

• все ограничения являются «мягкими»;

• часть ограничений является – «жесткими», другая часть – «мягкими».

Большинство реальных задач относится к третьему случаю. На пример, ограничения по безопасности химического производства от носятся к разделу «жестких», а ограничения на производительность и селективность часто могут быть отнесены к разделу «мягких».

Отметим, что при формулировании задачи интегрированного проектирования необходимым требованием является согласованность a, d и z в критериях для двух стадий «жизни» химико-технологической системы (требование реализуемости режимов функционирования ХТС).

Запишем постановку задачи ОЗИП в следующем виде (см. п. 2.1):

I = min M { C (a, d, z, )}, (2.50) a, d, z при связях в форме уравнений математической модели статики хими ко-технологической системы y = (a, d, z, ) (2.51) и ограничениях [ ] max max g j (a, d, z, ) 0, j J, j где g j (a, d, z, ) y j зад y j или [ ] max g j (a, d, z, ) y j зад y j 0, j J. (2.52) Заменим математическое ожидание с помощью квадратурной формулы некоторой суммой i C ( a, d, z, i ), M {C (a, d, z, ) iI i = 1 ;

I1 – множество аппрокси где i – весовые коэффициенты, iI мационных точек в области ;

y j, у j зад текущее и предельно допус тимое значения j-й выходной переменной.

Обозначим совокупность аппроксимационных точек i, i I через S1, а множество критических точек на -м шаге – через S2 ) = {l : j I 2l ) }. Алгоритм решения задачи (2.50), (2.51) основыва ( ( ется на методе внешней аппроксимации и его можно записать в сле дующем виде.

Алгоритм 2. Шаг 1. Полагаем число итераций = 1 и выбираем совокупность аппроксимационных точек S1, начальную совокупность критических точек S 20) и начальные приближения a ( 0), d ( 0), z ( 0).

( Шаг 2. Решаем вспомогательную задачу i C ( a, d, z, i ), min a,d, z iI g j (a, d, z, l ) 0, j = 1, m ;

l S 2 1) ;

l I 2 1), ( ( и определяем a ( ), d ( ), z ( ).

Шаг 3. Решаем m-задач max g j (a ( ), d ( ), z ( ), ), j = 1, m, и определяем m точек l ( ), l = 1, m.

Шаг 4. Образуем множество новых критических точек на -й итерации { } R ( ) = l ( ) : g j (a ( ), d ( ), z (), l ( ) ) 0.

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

Шаг 5. Формируем новое множество критических точек ( +1) = S 2 1) R ( ) и, полагая : = + 1, переходим к шагу 2.


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

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

Предполагаем на первом этапе, что функции gj выпуклы. В этом слу чае решение задачи max g j (a ( ), d ( ), z ( ), ), j = 1, m находится в одной из вершин параллелепипеда. В начальное множе ство критических точек S2(0) включается некоторое количество угло вых точек куба, а на шаге 3 рассчитываются значения функций g j (a ( ), d ( ), z ( ), ), j = 1, m во всех угловых точках куба, не принадлежащих множествам S2() и S1. Среди этих точек выбираются m точек, в которых функции g j (a ( ), d ( ), z ( ), l ), j, l = 1, m принима ют наибольшие значения. Далее формируем новое множество крити ческих точек S 2 +1) = S 2 ) R ( ) и переходим к шагу 2 алгоритма 2.2.

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

Имеются следующие трудности при решении задачи ДЗИП: 1) недиф ференцируемость функции h (a, d, ) по переменной и функций гиб кости 1 (a, d ) и 2 (a, d ) по переменной d;

2) многоэкстремальность задачи ДЗИП, связанная с невыпуклостью (невогнутостью) функций C (a, d, z, ), g j (a, d, z, ) и видом функции h (a, d, ).

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

Запишем дискретный вариант ДЗИП1:

i C ( a, d, z i, i ), С1= min a, d, z i iI g (a, d, z i, i ) 0, i I1, 1 (a, d ) = max min max g j (a, d, z, ) 0.

zZ jJ Рассмотрим метод внешней аппроксимации и метод разбиений и границ для решения дискретного варианта ДЗИП1 [31, 152]. Перепи шем эту задачу в виде i C ( a, d, z i, i ), C1 = min a, d, z i iI g (a, d, z i, i ) 0, i I1, h ( a, d, ) 0,, i и пусть [a, d, z ] решение этой задачи. Обозначим бесконечное множество точек, содержащихся в области, через S и множество точек p, которым соответствуют активные ограничения в точке ре шения задачи, через S A.P :

{ } S A.P = p : h(a, d, p ) = 0;

p. (2.53) i Будем называть эти точки активными. Решение [a, d, z ] задачи ДЗИП1 есть решение (локальный минимум) следующей задачи:

i C ( a, d, z i, i ), C1 = min a, d, z i Z iI g (a, d, z i, i ) 0, i I1, h(a, d, p ) 0, p S A.P.

Нижняя граница для ДЗИП1. Введем некоторое произвольное множество S 2 ) = { l : l I 2 ), l } точек из области неопределен ( ( ности, где I 2 ) – множество индексов точек в S 2 ) ( I1 I I 2 ) = ), и ( ( ( – номер итерации алгоритма решения ДЗИП1. Точки множества S2 ) будут называться критическими точками. Далее рассмотрим задачу ( i C ( a, d, z i, i ), C1L,( ) = min (2.54) a, d, z i iI g (a, d, z i, i ) 0, i I1, (2.55) h(a, d, l ) 0, l S 2 ).

( (2.56) C1L,( ) Эта задача имеет следующие свойства: 1) величина являет ся нижней границей оптимального значения целевой функции задачи ДЗИП1;

2) добавление точек к множеству критических точек не ухуд шает нижнюю границу;

3) если множество критических точек i, при надлежащих множеству S 2 ), достаточно плотно покрывает область, ( то решение задачи (2.54) – (2.56) достаточно близко к решению задачи ДЗИП1: С1 С1L,( ), где – достаточно малая положительная вели чина (точность решения ДЗИП1);

4) если решение a ( ), d ( ) задачи (2.54) – (2.56) удовлетворяет условию 1 (a, d ) 0, то a ( ), d ( ) есть решение дискретного варианта задачи ДЗИП1.

Задачу (2.54) – (2.56) перепишем в виде i C ( a, d, z i, i ), C1L,( ) = min (2.57) a, d, z i, z l iI g j (a, d, z i, i ) 0, j = 1,..., m, i I1, (2.58) g j (a, d, z l, l ) 0, j = 1,..., m, l S 2 ).

( (2.59) Сравнительный анализ задач (2.54) – (2.56) и (2.57) – (2.59) пока зывает, что прямое решение задачи (2.54) – (2.56) требует использова ния методов недифференцируемой оптимизации, поскольку функция h(a, d, ) недифференцируема. Кроме того, алгоритм решения этой задачи должен быть двухуровневым, так как в каждой точке [a, d, z i ] необходимо решать задачу нелинейного программирования 1 (a, d ) = max h(a, d, ), где h(a, d, ) = min max g j (a, d, z, ). Задача jJ z (2.57) – (2.59) имеет больше поисковых переменных, однако она явля ется более простой задачей: 1) она является одноуровневой задачей оптимизации;

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

Верхняя граница для ДЗИП1. Пусть область разбита на N i( ) ( i( ) = { : L,i,( ) U,i,( ) }), подобластей причем 1 ) = (N ), где номер итерации алгоритма решения зада U...U ( чи ДЗИП1. Заменим в задаче ДЗИП1 одно ограничение по гибкости N ограничениями 1,i (a, d ) 0 (i = 1,..., N ) [6]:

U i C ( a, d, z i, i ), С1,( ) = min U (2.60) a, d, z i iI g j (a, d, z i, i ) 0, j = 1,..., m, i I1, (2.61) 1,1 (a, d ) 0,..., 1, N (a, d ) 0, U U (2.62) где функция 1,i (a, d ) описывается формулой (2.29). В частности, ко U N = 1, ограничение 1 (a, d ) гда заменяется ограничением 1 (a, d) 0.

U Подставляя выражение (2.29) вместо 1,i (a, d ) в неравенства U (2.62) и используя ряд теорем из работ [31, 150], задачу (2.60) – (2.62) можно преобразовать к виду i C ( a, d, z i, i ), С1,( ) = min U (2.63) i l a, d, z, z iI g j (a, d, z i, i ) 0, j = 1,..., m, i I1, (2.64) max g j (a, d, z l, ) 0, j = 1,..., m, l = 1,..., N. (2.65) l ) ( Пусть [a ( ), z i,( ), z l,( ) ] решение задачи (2.63) – (2.65), тогда [a ( ), z i,( ) ] решение задачи (2.60) – (2.62).

j, (l ) Обозначим через решение задачи max g j (a ( ), d ( ), z l, ( ), ). Будем называть точку j, (l ) активной, ес () l ли соответствующее неравенство (2.65) является активным в точке решения задачи (2.63), т.е. g j (a ( ), d ( ), z l,( ), j,(l ) ) = 0.

Введем множество активных точек S AP задачи (2.63) – (2.65):

().

{ } S AP = l, j : g j (a ( ), d ( ), z l, ( ), j, (l ) ) = 0, l = 1,..., N, j = 1,..., m.

().

Будем называть область l( ) активной, если ей соответствует хо тя бы одно равенство g j (a ( ), d ( ), z l,( ), j, (l ) ) = 0. Активной области с номером l соответствует условие 1,l (a, d ) = 0. Ясно, что число ак U тивных областей не может быть больше размерности вектора конст руктивных параметров. Используя теорему П.6 [31, 150], можно запи сать задачу (2.63) – (2.65) в виде i C ( a, d, z i, i ), С1,( ) = min U a,d, zi, z l iI g j (a, d, z i, i ) 0, j = 1,..., m, i I1, g j (a, d, z l, j,(l ) ) 0, j = 1,..., m, l = 1,..., N, j, (l ) S AP.

().

Задача (2.60) – (2.62) имеет следующие свойства: 1) величина С1,( ) U является верхней границей оптимального значения целевой функции ДЗИП1;

2) дробление некоторых подобластей множества ( p ) не ухудшает верхнюю границу ДЗИП1, а в большинстве случаев даже ее улучшает;

3) если разбить область на достаточно малые подоблас ти i( ), то решение задачи (2.60) – (2.62) будет достаточно близко к C1,( ) = C1 ;

4) если размеры всех по U решению задачи ДЗИП1: lim r ( i ) ) ( добластей i( ) стремятся к нулю, то согласно свойству 3 задача (2.60) – (2.62) стремится к ДЗИП1 и ограничения (2.62) превращаются в ограни чения h(a, d, ) 0,. Поэтому множество активных точек S AP ().

задачи (2.60) – (2.62) стремится к множеству активных точек ДЗИП1.

Рассмотрим алгоритм вычисления верхней границы оптимального значения целевой функции ДЗИП1 (решения задачи (2.63) – (2.65)).

Обозначим через S 2k, p ) множество критических точек, принадлежа ( l щих подобласти l( ), где p номер итерации в алгоритме внешней аппроксимации решения задачи (2.63) – (2.65), номер итерации в алгоритме решения задачи ДЗИП1.

Множество S 2k, p ) будет образовываться автоматически при ре ( l шении задачи (2.60) – (2.62). Пусть S 2, p ) = S 2, p ) U S 2(2, p) U...U S 2(Np).

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

i C ( a, d, z i, i ), С1 L, ( ) = min U (2.66) a, d, z i, z l iI g j (a, d, z i, i ) 0, j = 1,..., m, i I1, (2.67) g (a, d, z l, l,q ) 0, l,q S 2l, p ), l = 1,..., N.

( (2.68), Задача (2.66) – (2.68) дает нижнюю границу оптимального значе ния целевой функции задачи (2.63) – (2.65). Сравним задачу (2.66) – (2.68) с задачей, дающей нижнюю границу оптимального значения целевой функции ДЗИП1 (см. задачу (2.57) – (2.59)). В задаче (2.66) – (2.68) один вектор z l соответствует всем точкам l,q, принадлежащим множеству S 2l, p ), в то время как в задаче (2.57) – (2.59) каждой точке (, i соответствует единственный вектор z i. Приведем описание алго ритма вычисления верхней границы оптимального значения целевой функции ДЗИП1 (решения задачи (2.63) – (2.65)).

Алгоритм 2. Шаг 1. Положим p = 1. Зададим множество подобластей l( ) (l = 1,..., N ) (число фиксировано, следовательно, число по добластей и их размеры остаются постоянными в этом алгоритме).

Зададим начальные значения z i,( 0), z l,( 0), a ( 0), d ( 0) (i I1, l = 1,..., N ) соответствующих переменных и достаточно малое число 0.

Шаг 2. Сформируем начальное множество S 2l,1) критических то (, чек для всех подобластей l (l = 1,..., N1 ) S2l,1) = {l, j,(1), j = 1,..., m : g j (a (0), d ( 0), z i,( 0), l, j,(1) ) 0}, (, где l, j,(1) есть решение задачи max g j (a ( 0), d ( 0), z l,( 0), ).

l ) ( Шаг 3. Решим задачу (2.66) – (2.68) и определим нижнюю грани цу С1,( ) величины С1,( ). Пусть [a ( p ), d ( p ), z i,( p ), z l,( p ) ] решение UL U этой задачи.


Шаг 4. Решим mN задач max g j (a ( p ), d ( p ), z l,( p ), ), где l ) ( { }. Обозначим через l( ) L,l, ( ) U,l, ( ) = : l, j,( p ) решение этой задачи при фиксированных [l, j ].

Шаг 5. Проверим выполнение условия g j (a ( p ), d ( p ), z l,( p ), l, j,( p ) ), j = 1,..., m;

l = 1,..., N.

Если все неравенства выполняются, решение задачи (2.63) – (2.65) получено, перейти к шагу 9;

в противном случае перейти к шагу 6.

Шаг 6. Создадим множества R l,( p ) (l = 1,..., N ) точек l, j,( p ), для которых ограничения нарушаются.

R l,( p) = {l, j,( p) : g j (a ( p), d ( p), z i,( p), l, j,( p) ) 0, j = 1,..., m, l = 1,..., N }.

Шаг 7. Сформируем новые множества критических точек S 2, p +l ) = S 2, p ) (l = 1,..., N ), U R l,( p) ( ( l l которые будут использоваться на следующей итерации.

Шаг 8. Положим p := p + 1 и перейдем к шагу 3.

Шаг 9. Сформируем и запомним множество активных точек S AP = {l, j,( p) : g j (a ( p), d ( p), z i,( p), l, j,( p) ) = 0, l = 1,..., N, j = 1,..., m }, ().

где p – есть номер последней итерации этого алгоритма.

Перейдем к описанию алгоритма решения ДЗИП1 в следующей постановке:

iC (a, d, z i, i ), С1 = min (2.69) a, d, z i iI g (a, d, z i, i ) 0, i I1, (2.70) max h(a, d, ) 0. (2.71) Задача (2.69) – (2.71) имеет вид задачи полубесконечного про граммирования и для ее решения можно использовать алгоритм внеш ней аппроксимации [40 – 42].

Вначале рассмотрим упрощенную задачу полубесконечного про граммирования:

С = min C ( x), (2.72) xX max g j ( x, ) 0, j = 1,..., m, { } где = : L U. Преобразуя ограничения задачи с помощью соотношения max g ( x) 0 g ( x) 0, x X, получим задачу xX С = min C ( x), xX g j ( x, ) 0, j = 1,..., m,, в которой число поисковых переменных конечно, а число ограничений бесконечно. Введем понятие множества критических точек S ( p ) = {i : i, i I ( p ) }, где I ( p ) множество индексов точек i.

Оно используется в алгоритме для построения некоторой аппроксима ции ограничений в задаче полубесконечного программирования (2.72).

Алгоритм вешней аппроксимации Шаг 1. Положим p = 1. Задаем начальное значение x ( 0) и на чальное множество критических точек S ( 0).

Шаг 2. Решаем задачу C ( p ) = min C ( x ), xX g j ( x, i ) 0, j = 1,..., m, i S ( p 1).

Пусть x ( p ) решение этой задачи.

Шаг 3. Решаем m задач max g j ( x ( p ), ) (2.73) и пусть (j p ) ( j = 1,..., m) решения этих задач.

Шаг 4. Проверяем выполнение условий g j ( x ( p ), (jp ) ) 0, j = 1,..., m. (2.74) Если все эти условия выполняются, то решение найдено.

Шаг 5. Сформируем множество R ( p ), которое содержит все точ ки j, в которых условия нарушаются:

( p) { } R ( p ) = (jp ) : g j ( x ( p ), (jp ) ) 0, j = 1,..., m.

Шаг 6. Сформируем новое множество критических точек S ( p ) = S ( p 1) U R ( p ).

Шаг 7. Полагаем p = p + 1 и переходим к шагу 2.

Этот алгоритм дает решение задачи (2.72), если операции на ша гах 2 и 3 выполняются в глобальном смысле.

На каждой итерации этот алгоритм выполняет две основные опе рации. Первая операция связана с определением нижней границы це левой функции задачи (2.72). Вторая операция связана с проверкой, является ли точка x ( p ) решением задачи (2.72). Здесь алгоритм решает задачу (2.73) m раз. Если условие (2.74) удовлетворяется, то решение задачи (2.72) найдено. В противном случае точки (jp ), в которых ус ловия (2.74) нарушаются, добавляются в множество S ( p ).

Заметим, что ОЗИП (2.50) с жесткими ограничениями имеет вид (2.72), поэтому для ее решения также можно использовать описанный алгоритм внешней аппроксимации. Существенный недостаток описан ного алгоритма состоит в необходимости решения m задач (2.73) на каждой итерации. Авторы работы [43] предложили заменить m задач максимизации (2.73) одной задачей максимизации функции F (g ) :

max F ( g ), 1 m g j F (g) = ln e, j =1 где – некоторый параметр. В работе доказано, что F ( x, ) max ( g j ( x)), 0, j lim F ( x, ) max ( g j ( x)).

j Теперь вернемся к решению задачи (2.69) – (2.71). Она имеет вид задачи полубесконечного программирования и для ее решения можно использовать описанный выше алгоритм внешней аппроксимации.

Поскольку max h(a, d, ) = 1 (a, d ), то на каждой итерации будет не обходимо вычислять функцию гибкости 1 (a, d ) (шаг 3 алгоритма внешней аппроксимации).

Совокупность аппроксимационных точек i, i I1, будем обозна чать через S1, а множество критических точек на -м шаге – через S 2 ) = { i : i I 2 ) }.

( ( Алгоритм 2. Шаг 1. Положим = 1. Задаем множество аппроксимационных точек S1, начальное множество критических точек S 20), начальные ( приближения a ( 0), d ( 0) и достаточно малое число 0.

Шаг 2. Решаем задачу (2.57) – (2.59) для определения нижней границы С1L,( ) оптимального значения величины C1. Пусть a ( ), d ( ) решение этой задачи.

Шаг 3. Находим значение 1 (a ( ), d ( ) ), и пусть ( ) – решение этой задачи.

Шаг 4. Если выполняется условие 1 (a ( ), d ( ) ) 0, то решение найдено. В противном случае находим точку ( ), для которой нару шаются ограничения g j (a ( ), d ( ), ( ) ) 0, j = 1,..., m, и переходим к шагу 5.

Шаг 5. Образуем новое множество критических точек S 2 +1), до ( U {() }.

бавляя точку ( ) к множеству S 2 ) : S 2 +1) = S 2 ) ( ( ( Шаг 6. Полагаем = + 1 и переходим к шагу 2.

Для вычисления значения 1 (a ( ), d ( ) ) функции гибкости можно использовать метод ветвей и границ. В статье [44] было предложено использовать на шаге 3 следующее приближение ~1 (a ( ), d ( ) ) значе ния 1 (a ( ), d ( ) ) :

1,( ), если 1,( ) 0 ;

U U U,( ) + 1,( ) L 1 (a ( ), d ( ) ) = 1 = 1 если 1,( ) 1,( ) ;

~ U L, 1,( ), если 1,( ), L L где 1,( ), 1,( ) являются верхней и нижней границей функции гиб L U кости 1 (a ( ), d ( ) ). В случае использования метода ветвей и границ величина 1,( ) определяется по формуле (2.35), в которой на каждой U итерации величины 1,i,( ) известны для каждой подобласти. Величина U 1,( ) получается с помощью формулы (2.42). Заметим, что если усло L вие 1,( ) U выполняется, то тем более выполняется условие ( ) ( ) 1 (a ) 0. Следовательно, если выполняется либо условие,d 1,( ) 0, либо условие 1,( ) 1,( ), где мало, то решение U L U ДЗИП1 получено. Выполнение условия 1,( ) (a, d ) 0 означает, что L имеет место условие 1 (a, d ) 0.

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

Далее рассмотрим метод разбиений и границ для решения ДЗИП1 [6, 7], представляющий двухуровневую итерационную проце дуру (рис. 2.3).

Верхний уровень i(v ) i(v ) С1, (v ), С1L, (v ) U Вычисления верхней С1, (v ) Нижний U и нижней уровень С1L, (v ) границ Рис. 2.3. Блок-схема метода разбиений и границ Он состоит из следующих блоков: 1) вычисление верхней грани цы С1,( ) (см. задачу (2.60) – (2.62)) целевой функции задачи U i C ( a, d, z i, i ), С1= min (2.75) a, d, z i iI g (a, d, z i, i ) 0, i I1, (2.76) 1 (a, d ) = max min max g j (a, d, z, ) 0 ;

(2.77) jJ z 2) алгоритм вычисления нижней границы С1L,( ) (см. задачу (2.57) – (2.59)) целевой функции задачи (2.75) – (2.77);

3) разбиение области неопределенности (правило разбиения);

4) формирование множества S 2 ) критических точек (правило выбора) в алгоритме вы ( числения нижней границы.

Верхний уровень алгоритма разбиения и границ (см. рис. 2.3) служит для проверки окончания процедуры решения и разбиения об ласти на подобласти. На нижнем уровне вычисляются верхняя С1,( ) и нижняя С1L,( ) границы.

U Сформулируем правила разбиения области на подобласти () i (i = 1,..., N ). Прямой путь состоит в систематическом дроблении всех областей на каждой итерации, пока все подобласти не станут дос таточно малыми (свойство 3 задачи (2.60) – (2.62)). Однако такая стра тегия разбиения неэффективна, так как в этом случае размерность за дач (2.57) – (2.59) и (2.63) – (2.65) достаточно быстро может стать очень высокой.

Рассмотрим более эффективный способ дробления, основанный на следующем эвристическом правиле: на -й итерации дробиться будут только те из областей i( ) (i = 1,..., N ), для которых соот ветствующие ограничения (2.61) задачи (2.60) – (2.62) будут активны.

Так как задача (2.63) – (2.65) эквивалентна задаче (2.60) – (2.62), то это эвристическое правило может быть сформулировано следующим обра зом: на -й итерации подобласть l( ) дробится, если выполняется условие j J max g j (a ( ), d ( ), z l,( ), ) = 0, где J = (1,..., m). Здесь l () ( ) l,( ) ] решение задачи (2.63) – (2.65).

[a,d,z В работе [31] показано, что дробление подобластей с неактивны ми ограничениями бесполезно, а удаление активных ограничений в задаче оптимизации улучшает оптимальное значение целевой функции.

В качестве множества S 2 ) можно использовать все критические ( точки, полученные при решении задачи (2.63) – (2.65) вычисления верхней границы целевой функции ДЗИП1. Однако в этом случае раз мерность задачи (2.57) – (2.59) может стать очень высокой. Необходи мо выбирать это множество таким образом, чтобы нижняя граница C1L (a, d ) была как можно ближе к оптимальному значению C1 целе вой функции ДЗИП1. Как было показано ранее, C1, ( ) (a, d ) стремится U к C1 при max r ( l( ) ) 0. Поэтому целесообразно, чтобы при l N задача (2.57) – (2.59) вычисления нижней границы стреми лась к задаче вычисления верхней границы:

i C ( a, d, z i, i ), С1,( ) = min U a, d, z i, z l iI g j (a, d, z i, i ) 0, j = 1,..., m, i I1, g j (a, d, z l, jl ) 0, j = 1,..., m, l = 1,..., N, jl S AP.

().

Следовательно, целесообразно выбрать множество S 2 ) следую ( щим образом: S 2 ) = S AP. При этом можно уменьшить число точек, ( ().

используемых для формирования S 2 ). Действительно, при N ( все точки одной подобласти стягиваются в одну точку. Следовательно, во множество S 2 ) можно включать только одну активную точку из ( всех точек, принадлежащих подобласти l( ). В связи с этим образуем ~( ) множество S AP из множества S AP следующим образом. В каждой ()..

l подобласти l( ) из всех активных точек S A,.(P ) будем оставлять только ~ ~ одну, которую будем обозначать через l. Кроме того, в качестве l можно выбрать среднее арифметическое из всех активных точек по добласти l( ) :

~l l, j, = l S A,.(P ) jJ l, ( ) A где S A,.(P число точек во множестве S A,.(P ) ;

J A,( ) множество номе l ) l l ров активных точек подобласти l( ).

~( ) Множество S AP будет состоять из образованных таким образом.

~l точек. В качестве множества S 2 ) будем использовать множество ( ~( ) ~( ) S AP, т.е. S 2) = S AP.

(..

Рассмотрим метод разбиений и границ решения задачи ДЗИП1.

Введем множество L( ) подобластей i( ) следующим образом:

{ } L1 ) = i( ), (i = 1,..., N ) : r (i( ) ) 1, ( где 1 заранее заданное число.

Алгоритм 2. Шаг 1. Положим = 1. Зададим начальное множество подобла l(1) (l = 1,..., N ), множество стей аппроксимационных точек { } S1 = i : i I1, начальное множество критических точек S 20), началь ( z i,( 0), z l,( 0), a ( 0), d (0) (i I1, l = 1,..., N1 ) ные значения соответствую 1 щих переменных, число и достаточно малые числа 1 0, 2 0, 2 0, (2 1, 1 2 ). C1,(0) = a, U Положим C1L,( 0) = a, где a достаточно большое число (a C1 ).

Шаг 2. Вычислим верхнюю границу С1,( ) решением задачи U с помощью алгоритма 2.3. Пусть (2.60) – (2.62) [a ( ), d ( ), z i,( ), z l,( ) ] (i I1 ), (l = 1,..., N ) решение этой задачи.

{ } Шаг 3. Определим множество Q ( ) = l( ) : l I Q ) подобластей ( i( ), которым соответствуют активные ограничения в задаче (2.60) – (2.62):

1,l (a ( ), d ( ) ) = 0, l( ) Q ( ).

U Шаг 4. Если Q ( ) пустое множество, то решение ДЗИП1 найде но, алгоритм заканчивает свою работу.

Шаг 5. Если выполняется условие С1,( 1) С1,( ) 1 С1,( ), U U U то решение ДЗИП 1 найдено, алгоритм заканчивает свою работу.

В противном случае проверяем выполнение условия С1,( 1) С1,( ) 2 С1,( ), U U U и если оно нарушается, то переходим к шагу 8, в противном случае переходим к шагу 6.

Шаг 6. Находим нижнюю границу С1L,( ), решая задачу (2.57) – (2.59). Здесь мы полагаем S 2 ) = S AP (см. шаг 9 алгоритма 2.3).

( ().

Шаг 7. Если выполняется условие C1,( 1) C1,( ) 1 C1,( ), U U U то решение ДЗИП1 найдено, алгоритм заканчивает свою работу.

Шаг 8. Если выполняется условие r ( i( ) ) 2, i = 1,..., N, то решение ДЗИП1 найдено, алгоритм заканчивает свою работу.

Шаг 9. Если выполняется условие r (l( ) ) 1, l = 1,..., N, то переходим к шагу 11.

Шаг 10. Сформировать множество L( ). Найти множество V ( ) подобластей i( ), принадлежащих одновременно множествам L( ) и Q ( ), т.е.

V ( ) = L( ) I Q ( ).

Разбить каждую подобласть l( ) V ( ) на две подобласти таким образом, что l( +1) U l(+1) = l(). Образовать новое множество 1 ( +1) подобластей из старого множества ( ), заменяя подобласть l( ) новыми подобластями l( +1) и l( +1). Переходим к шагу 12.

1 Шаг 11. Положим 1 = 1 / 2 и переходим к шагу 9.

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

Алгоритм 2.5 позволяет получить локальный минимум.

Определение начального разбиения области на подобласти на шаге 1 является непростой задачей. Конечно, начальное разбиение может состоять из всей области. Однако в этом случае задача (2.60) – (2.62) может не иметь решения. В связи с этим на шаге 1 мож но вычислять значение индекса гибкости 1 ХТС. Эта величина нахо дится решением задачи 1 = min max min max g j (a, d, z, ) jJ a, d z или 1 = min u, (2.78) a A, d D 1 (a, d ) = max min max g j (a, d, z, ) u. (2.79) jJ z Ясно, что если 1 0, то ДЗИП1 не имеет решения. Если 1 0, то, решая задачу (2.78), (2.79), мы получаем некоторые значения a, d и некоторое разбиение области, для которой задача ДЗИП i C ( a, d, z i, i ), С1= min (2.80) a, d, z i iI g (a, d, z i, i ) 0, i I1, (2.81) 1 (a, d ) = max min max g j (a, d, z, ) 0 (2.82) jJ z имеет решение. Множество подобластей, найденное решением задачи (2.78),(2.79), может быть использовано как начальное разбиение об ласти на подобласти i(1) на шаге 1 алгоритма 2.5. Задача (2.78), (2.79) представляет собой частный случай ДЗИП1. Поэтому для ее ре шения можно использовать алгоритм 2.5.

До сих пор мы предполагали, что аппроксимационные точки в дискретном варианте ДЗИП1 (см. задачу (2.80)–(2.82)) задаются поль зователем, и что число этих точек невелико. Здесь мы рассматриваем случай, когда необходимо решать задачу (2.14) – (2.15). В данном слу чае главной проблемой является вычисление значения критерия, т.е.

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

Первый подход основывается на использовании квадратурных формул Гаусса [35]. Предположим для простоты изложения, что век тор имеет только две компоненты (т.е. область прямоугольник);

тогда, если все неопределенные параметры независимы, мы должны вычислить интеграл U U C (a, d, z(), () P())d1d 2.

I= (2.83) LL 2 Введем p точек (1 ) и ( k ) на интервалах [1, 1 ] и [2, U ], k LU L 2 соответственно:

1 = 1 + k1, k = 0, 1,..., ( p 1), k L k = 2 + k 2, L (1 1 ) ( U 2 ) U L L где 1 =, 2 = 2.

p p Пересечения прямых (рис. 2.4): 1 = 1, i = 0, 1,..., ( p 1) i и 2 = 2, j = 0, 1,..., ( p 1), образуют p 2 точек пересечения (узлов) j ij = (1, 2 ) (i = 0,..., p 1, j = 0,..., p 1), 2 = 2 + j 2, i j j L где 1j = 1 + j1, и квадратурная формула Гаусса приближенного значе L ния интеграла (2.83) имеет вид (1 1 ) (U 2 ) p 1 i p 1 j U L L I 1 2 C (a, d, z ij, 1, 2 ), j i 2 2 i =0 j = где z ij управляющая переменная, соответствующая точке ij. Квад ратурная формула легко обобщается на п-мерный случай. В этом случае для аппроксимации многомерного интеграла в задаче (2.14) – (2.15) n потребуется p узловых точек, где n размерность вектора.

Заменяя в задаче (2.14) – (2.15) целевую функцию ее гауссовым при ближением, получим дискретный вариант ДЗИП1 (2.80) – (2.82). Здесь узловые точки используются в качестве аппроксимационных точек.

{ } Следовательно, множество S1 = i : i I1 аппроксимационных точек n в задаче (2.80)–(2.81) будет состоять из p точек. Для вычисления i +1 j ij v j L 1 L i v Рис. 2.4. Узловые точки двумерного интеграла целевой функции в задаче (2.80)–(2.81) мы должны вычислить значе n ние функции С (a, d, z, ) p раз. Размерность задачи (2.57) – (2.59) n вычисления нижней границы будет равна nd + ( p + N c. p )n z, где nd и n z размерности векторов d и z, соответственно, N c. p число кри тических точек. Очевидно, что решение практических задач большой размерности может потребовать огромных вычислительных затрат.

Второй подход к вычислению интегралов основывается на проце дуре Монте-Карло и близких к ней процедурах (латинского гиперкуба и последовательности проб Хаммерслея (HSS) [45, 46]. Для этих мето дов известно, что число аппроксимационных точек, необходимых для вычисления интеграла с заданной точностью, не зависит от размерно сти вектора и техника HSS наиболее эффективна по сравнению со всеми подобными подходами. При этом даже метод HSS требует не скольких сотен аппроксимационных точек для получения разумной точности вычисления интеграла.

Техника Монте-Карло для оценки интеграла в задаче (2.14) – (2.15) достаточно проста. Пусть имеется последовательность, состоя щая из N векторов i, имеющих плотность распределения вероятно стей P(). Тогда интеграл можно вычислить по формуле C (a, d, z (), P())d N C (a, d, z( ), ).

i i iI В работе [32] в качестве аппроксимационных точек использовали номинальную точку и все критические точки, получаемые в итерациях этого метода. Весовой коэффициент, соответствующий номинальной точке, выбирался равным 0,5, остальные коэффициенты выбирались равными 0,5 / (n – 1).

Видно, что выбор множества аппроксимационных точек в данном случае достаточно субъективен. В связи с этим рассмотрим более объек тивный подход к построению целевой функции в ДЗИП1. В каждой за даче можно выделить три характерные точки. Первая точка – номиналь ная точка N. Если распределение (неизвестное) является симметрич ным, то это наиболее вероятная точка. Вторая и третья точки соответст вуют точкам с наилучшим и наихудшим значениями функции при фик сированных a, d, z. На основе этих точек можно сформулировать три целевые функции и три внутренние задачи оптимизации ДЗИП:

C1 (a, d ) = min C (a, d, z, N ), z g j (a, d, z, ) 0, j = 1,..., m.

C2 (a, d ) = min min C (a, d, z, ), z g j (a, d, z, ) 0, j = 1,..., m.

C3 (a, d ) = min max C (a, d, z, ), z g j (a, d, z, ) 0, j = 1,..., m.

На основе этих внутренних задач сформируем следующие три ва рианта ДЗИП:

Ci = min Ci (a, d ), i = 1, 2, 3, a,d 1 (a, d ) 0.

Существуют неравенства: C2 C C3, где C получается решением ДЗИП С1= min i C ( a, d, z i, i ), a,d, z i Z iI g (a, d, z i, i ) 0, i I1, 1 ( a, d ) = max min max g j ( a, d, z, ) zZ jJ для любого набора аппроксимационных точек и весовых коэффициен тов. Таким образом, решив задачу ДЗИП для вариантов i = 2, 3, можно оценить диапазон, в котором будет лежать оптимальное значение це левой функции ДЗИП1.



Pages:     | 1 | 2 || 4 | 5 |   ...   | 8 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.