CC BY
121
№ 3 (45) 2013
А. А. Емельянов, докт. экон. наук, профессор Национального исследовательского университета «МЭИ», г. Москва
Планирование экстремальных экспериментов с имитационными моделями
В случае если набор выходных параметров модели, получаемых автоматически, не устраивает экспериментатора, причем поверхность отклика, где можно найти экстремальные значения, неизвестна, и при этом необходимо оптимизировать функционирование поведенческой модели реального процесса путем подбора входных параметров-факторов, нужно использовать специальные методы проведения научного эксперимента [10,11].
введение
Имитационная модель, независимо от выбранной системы моделирования Actor Pilgrim, AnyLogic, GPSS World или ReThink (см. тематический обзор [4]), позволяет автоматически получить численные значения двух первых моментов любой получаемой в процессе моделирования величины и информацию о ее распределении. Кроме того, для любого входного параметра-фактора, интересующего экспериментатора, можно задавать обоснованно выбранное им распределение либо передавать распределение, получаемое в исследуемой или иной модели.
При моделировании сложных экономических процессов разработчики моделей исходят из того, что за время прогона реализуются от сотен тысяч до миллиарда испытаний поведенческих процессов, и поэтому первые моменты искомых параметров получаются с должной достоверностью. Однако при моделировании рисковых процессов, когда из миллионов поведенческих процессов нужно получить хотя бы несколько таких, которые привели к рисковой ситуации, возникают следующие существенные трудности.
1. Астрономическое время одного прогона модели на современном компьютере может продолжаться много часов, а иногда превысить сутки, по нижеприведенным причинам.
2. Во время многодневного (как правило) эксперимента по исследованию рисковых поведенческих процессов необходимо провести имитационное расследование цепочек причин, приводящих к рисковым событиям. Однако в ходе моделирования удается наблюдать сравнительно немного таких процессов, поскольку далеко не при каждом опыте возникает модельное рисковое событие. Соответственно, и выводы, как правило, делаются на основании изучения в лучшем случае десяти таких процессов. Эта трудность является следствием предшествующей.
3. Датчики случайных величин в современных моделирующих пакетах оставляют желать лучшего, так как их период псевдослучайных последовательностей обычно не превышает 0,5 млрд. Поэтому, если рисковая ситуация не обнаружена по истечении такой последовательности, то простое увеличение числа испытаний в несколько раз все равно не приведет к получению результата, так как моделирование будет повторяться на одном наборе псевдослучайных чисел. Исключением являются пакеты MatLab, у которых длина периода = 4,72236197х1021, и Actor Pilgrim, где длина периода = 1,62259276х1032 [6]. Однако у двух вышеназванных пакетов назначение различное:
— MatLab предназначен для моделирования главным образом физических, техни-
№ 3 (45) 2013
ческих и технологических процессов, а также рисков, связанных с технической и технологической надежностью и неисправностями оборудования;
— Actor Pilgrim можно использовать для тех же целей, что и MatLab. Но этот пакет имеет следующие отличительные функциональные возможности:
1) создание моделей сложных поведенческих процессов в экономике, с учетом финансовой динамики и движения денежных средств;
2) моделирование рисковых поведенческих процессов, приводящих к банкротству, к инвестиционным неудачам;
3) симуляция клиринговых процессов в банковской сфере;
4) имитация территориально-региональных экономических, народно-хозяйственных и экологических процессов с использованием топографической информации, в том числе процессов гражданской защиты.
Целью данной статьи является демонстрация тех возможностей проведения имитационных экспериментов с обоснованием уровня доверия к полученным результатам, которые получаются с помощью моделей, созданных в Actor Pilgrim [5].
способы оценки первых двух моментов выходных параметров при моделировании сложных процессов без рисковых событий
Определим базовые термины1. Назовем испытанием один прогон модели, позволяющий получить конкретный численный результат.
Термином опыт будем называть серию прогонов модели, которая дает возможность получить как минимум два первых момента интересующих нас выходных параметров.
1 Здесь и далее в статье используется терминология, которая применялась некоторыми авторами при обсуждении вопросов планирования эксперимента
[7, 10, 11].
Экспериментом назовем серию целе- | направленных опытов, позволяющих путем г§ варьирования параметров-факторов найти ¿g такое сочетание этих параметров, которое ^ с заданной достоверностью позволяет по- ^ лучить искомое значение выходного параметра отклика.
Рассмотрим получение первого и второго момента произвольного интервала времени. Эта величина может измеряться двумя способами. Далее полагаем, что имеем дело с имитационными моделями, созданными в системе Actor Pilgrim.
Первый способ. Искомая величина является интервалом задержки актора в узле. Дисперсия получается как произведе-
2
ние квадрата математического ожидания m на коэффициент вариации в квадрате с2, который также подсчитывается автоматически в этом узле.
Второй способ. Искомая величина не является параметром задержки в узле модели. Однако можно применить два методических приема и связать ее следующими параметрами:
1) интервал пребывания клапана key, дополнительно введенного в модель, в запертом состоянии (математическое ожидание этого интервала и квадрат коэффициента вариации вычисляются в узле key автоматически);
2) время жизни дополнительно сгенерированного транзакта в узле create, помещенного в запертый клапан key, который в нужный момент открывается и направляется в дополнительный терминатор term (математическое ожидание и квадрат коэффициента вариации времени жизни определяется в узле term автоматически).
Для получения вида неизвестного закона распределения, если нет оснований для предложения и доказательства гипотезы об этом законе, существует простейший прием. Интересующий интервал возможных значений переменой x, которая имеет произвольный смысл (денежная сумма, объем партии товара и др.), делится на k равных интервалов:
№ 3 (45) 2013
I
!
1
it
0
и
is
1
to
I
та
i
*
I
s
СО
0
1 ¿1
(x0, Xi], (Xi, x2], ..., (xk_ 1, xk].
В модели объявляется массив переменных с фиксированной точкой: long p [k].
Во время прогона модели частоты появления значений x в этих интервалах значений подсчитываются в соответствующих элементах массива p [k]. Вид закона определяется в виде ступенчатой функции. Поэтому получение доверительного интервала значений измеряемой величины x или проверка гипотезы о равенстве математического ожидания M(x) заданному значению d = const не вызывает затруднений.
По поводу достоверности результатов можно отметить следующее: при моделировании сложных процессов без рисковых событий, когда устанавливается стационарный режим в рамках одного испытания, вопрос о доверительных интервалах не возникает, так как число испытаний обычно составляет от 106 до 108 в пределах одного опыта. В результате, как правило, математические ожидания и среднеквадратичные отклонения имеют отклонения в четвертом знаке. Однако в случае сомнений число испытаний увеличивается на порядок.
оценки достоверности эксперимента по выборке поведенческих сценариев с рисковыми событиями из генеральной совокупности
Итак, в выборку необходимо включать результаты только таких опытов, в которых замечено хотя бы одно рисковое событие. Зачастую в них предстоит выполнить не менее 1010^1011 испытаний, причем последовательности в серии опытов не совпадают (они всегда новые). В моделях для этого используются генераторы сверхдлинных последовательностей случайных величин. Поскольку в основном будут иметь место опыты, где рисковых событий не происходит, объем выборки п будет существенно меньше 100: выборка объемом 10 опытов — хороший результат.
Доверительные интервалы. Промежуток, которому принадлежит оцениваемый параметр 0 с вероятностью > q, называют доверительным интервалом, q — коэффициентом доверия, а величину, равную 1- q, — уровнем значимости [3]. Введем следующие обозначения: вп — несмещенная оценка параметра 0 по выборке объемом п, е — допустимое максимальное отклонение. Тогда о справедливости условия
Р (|0-0ч|<е)> q,
означающего, что 0 е ( -е, ©п -е) с вероятностью > q, можно судить с помощью неравенства Чебышева [12], но это даст только грубую оценку.
Однако соответствующий «рецепт» очевиден. Если 0 — математическое ожидание X, а вп — его несмещенная оценка, то
Р (М„|<е)> 1-^.
Практический способ действия на этой основе заключается в следующем. Задается коэффициент доверия
q = 1-
D(0n
откуда е=
D(Qn)
1-q
что определяет доверительный интервал
(-е е„-е).
Здесь в некотором роде заложено противоречие, поскольку на практике обычно имеется реализация выборки и более — ничего. Поэтому в получаемых неравенствах неизвестное оценивается через неизвестное. Дисперсию 0(еп) приходится определять по той же выборке. Однако противоречие снимается, если оценки состоятельны. Тогда 0(еп) определяется с небольшой ошибкой А, и
е =
D0J 1 - q
+ o(q),
где о^) << q — малая функциональная величина, т. е. влиянием ошибки при определении дисперсии можно пренебречь.
78
2
е
№ 3 (45) 2013
Если речь идет о достаточно длинных выборках, то можно опираться на предельные теоремы о нормальности распределения ошибок при усреднении, что дает более точные оценки.
Пусть, например, оценивается вероятность р некоторого события А, где Хк принимает значения ноль/один при к-м испытании в выборке Х1, Х2, ..., Хп: Хк = 1 означает наступление события А. Если для оценки используется среднее
р = Х1 + Х2 +... + Хп
п
то, очевидно,
М (Рп) = Р, О (Рп) = ^.
При больших выборках п, в силу предельных теорем
Р - Pj<ei
Рк (1-Рп)
п
= 2Ф(е),
где Ф(е) — функция распределения нормального закона распределения, откуда получается необходимая связь между крайними точками доверительного интервала и уровнем значимости.
Очевидно, что строгое решение последнего уравнения занимает много места, но оно приходится кстати, например в случае, когда исследователю-теоретику требуется представить масштабный отчет о проделанной работе. Для практических имитационных экспериментов доверительный интервал приближенно равен
(Рп - е°' Рп+ео)'
где о =
fiK).
Все это хорошо и точно работает, когда выборка достаточно велика (практически П = 102). Однако для малых п приходится начинать эксперимент с предположениями о распределении Я. Бернулли, что в отсутствии возможности воспользоваться фор-
мулой Стирлинга [12] приводит к дополнительным компьютерным вычислениям, но, по мнению автора данной статьи, не является недостатком [3].
введение в теорию планирования экстремального модельного эксперимента
Различают два вида эксперимента: пассивный и активный.
Пассивный эксперимент обычно выполняется на базе действующего экономического, технологического или иного процесса, который не допускает активного вмешательства экспериментатора. Этот метод малозатратный, но требует большого времени. При пассивном эксперименте исследователь только ведет наблюдение за процессом, за изменением его входных и выходных параметров. По результатам наблюдений затем делается вывод о влиянии входных параметров на выходные. Пассивный эксперимент практически невозможен для исследования рисковых процессов, связанных с экономическими реформами, банкротством, гражданской защитой или военными действиями.
Активный эксперимент проводится главным образом в лабораторных условиях, где экспериментатор имеет возможность изменять входные характеристики по заранее намеченному плану. Такой эксперимент быстрее приводит к цели и именно к нему применимы идеи планирования экстремального эксперимента.
Планирование имитационного эксперимента — это сложная управленческая задача поиска и определения области, в которой находится оптимальная в каком-то смысле, с точки зрения экспериментатора, точка.
Пла ниро в ан и е экс пе р и ме нт а можно рассматривать как кибернетический подход к организации и проведению экспериментальных исследований сложных объектов и процессов. Основная идея метода состоит в возможности оптимального
0
1
79
управления экспериментом в условиях неопределенности, что родственно тем предпосылкам, на которых базируется кибернетика. Целью большинства исследовательских работ является определение оптимальных параметров сложной системы или оптимальных условий протекания процесса, например:
— определение параметров инвестиционного проекта в условиях неопределенности и риска;
— выбор конструкционных и электрических параметров энергетической установки, обеспечивающих наиболее выгодный режим ее работы;
— получение максимально возможного выхода продукта с заданными параметрами из химического реактора путем варьирования температуры, давления, соотношения реагентов и используемых катализаторов.
При решении задач такого рода приходится учитывать влияние большого количества факторов, часть из которых не поддается регулированию и контролю, что чрезвычайно затрудняет полное теоретическое ^ исследование задачи. Поэтому идут по пу-| ти установления основных закономерно-Ц стей с помощью проведения серии экспе-| риментов. Методы эмпирического поиска | оптимального решения долгое время оста-| вались неформализованными. Исследова-§ тель выбирал ту или иную схему постановки I эксперимента (стратегию), базируясь только 15 на своем опыте и интуиции. Однако начиная § со второй половины XX в. стала усиленно ^ развиваться математическая теория экстремальных экспериментов, которая помога-| ет экспериментатору выбрать оптимальную « стратегию [11]. Основными показателями Ц оптимальности при этом является уменьше-| ние числа экспериментов при обеспечении Ц той же точности результатов исследования [з или сохранение числа экспериментов при увеличении точности. Существенные упро-| щения при этом достигнуты в методах об-§ работки результатов эксперимента. Иссле-| дователь получил возможность путем неё сложных вычислений выражать результа-
ты эксперимента в удобной для их анализа и использования форме.
В общем случае объект исследования можно представить как некоторый «черный ящик» (рис. 1), на входе которого действуют управляющие параметры хп, I = 1, 2, ..., N, и параметры неконтролируемых возмущений г,, I = 1,2, ..., L. Выходом объекта исследования являются показатели эффективности, качества или какие-либо другие характеристики объекта ут, I = 1,2, ..., М.
Рис. 1. Объект исследования в виде «черного ящика»
Например, в задачах анализа бизнес-процессов под переменными хп можно понимать номенклатуру материалов для выпуска изделий и финансирование для запуска производства принципиально новых товаров и услуг вместо тех, реализация которых стала затруднительной, под помехами г, — спрос на рынке во времени, ставки налогообложения, а под выходом ут — финансово-хозяйственные показатели.
Если рассмотреть зависимость одной из характеристик системы ут (хп) как функцию только одной переменной хп (рис. 2), то при фиксированном значении хп будем получать различные значения ут (хп). Разброс значений ут в данном случае определяется как ошибками измерения хп, так и влиянием помех г,. График зависимости ут (хп) определяет корреляционную связь величин ут и хп, которая может быть получена по результатам эксперимента с помощью методов математической статистики. Вычисление таких зависимостей при боль-
№ 3 (45) 2013
шом числе входных параметров хп и существенном влиянии помех г, и является основной задачей исследователя-экспериментатора. Причем чем сложнее задача, тем эффективнее становится применение методов планирования эксперимента.
Математически задача установления взаимосвязей оптимизируемого процесса формулируется следующим образом: требуется получить некоторое представление о неизвестной функции
Y = 1 (X),
где Y = (у; у2; ...; ук) — вектор выходных параметров процесса, подлежащий оптимизации;
X = (х1; x2
Далее для простоты предположим, что Y = (у1) — скаляр, т. е. К = 1 (заметим, что в случаях, когда К > 1, всегда можно ввести скалярное отображение вектора Y на некую координатную ось — например, модуль).
Естественно, и в этом случае аналитическое выражение функции отклика неизвестно. Наиболее удобным оказалось представление функции отклика в виде полинома, поскольку разложение функции в степенной ряд эквивалентно представлению ее рядом Тейлора [11]:
У = b 0 + %brXn +
n, u xrxu + Ibr
, x n+.
(1)
xN) — вектор входных пара-
мые переменные, которые можно варьировать при постановке экспериментов, в связи с чем они названы факторами.
Рис. 2. Пример усреднения результатов эксперимента
Конкретное значение точки с координатами Yk называется откликом, а Y = 1(Х) — функцией отклика. Геометрический образ соответствующей функции отклика Y = 1( X) называется поверхностью отклика. Пространство с координатами х1, х2, ..., хы называется факторным пространством.
где у — приближение Y = (у1), называемое уравнением регрессии, получаемое с помощью этого полинома, поскольку математический вид функции 1(Х) неизвестен; п, и = 1, 2, ..., N; Ь0 , Ьп , Ьпи , Ьпп , ..., — коэффициенты регрессии, т. е. приближения коэффициентов ряда Тейлора, получаемые только в результате эксперимента2 с хорошим уровнем достоверности, поскольку эти коэффициенты аналитически найти невозможно, так как они определяются через частные производные неизвестной функции отклика 1 (X).
Полином вида (1) называется уравнением регрессии.
Для определения численных значений выборочных коэффициентов регрессии используется так называемый регрессионный анализ, основанный на методе наименьших квадратов. Обычно в регрессионном анализе полагается, что выполняется ряд предпосылок. Главными считаются следующие предпосылки.
1. Результаты наблюдений у1, у2, ..., ум — независимые, нормально распределенные случайные величины. Если эта предпосылка не удовлетворяется, то коэффициенты ре-
2 Получение коэффициентов регрессии в модельном эксперименте достаточно подробно показано в [7, 11].
81
0
1 ¡1
r=1
n=1
метров, причем x1, x2, ..., xN — независи-
№ 3 (45) 2013
грессии наити можно, однако ничего нельзя будет сказать об эффективности метода, т. е. нельзя оценить точность уравнения регрессии. Если ут не подчиняется нормальному распределению, то стараются подобрать такую функцию преобразования, чтобы переИти от ут к новоИ случайной величине qm = ф( ут), распределенной «почти» нормально. Например, для многих асимметричных распределений делается замена нормального распределения на логнормаль-ное [6]: qm = 1п( Ут).
2. Если производить многократные и повторные наблюдения над величиной ут на некотором определенном наборе значений, то получим дисперсию D(ут) = а^у , которая не будет зависеть от математического ожидания М( ут), т. е. не будет отличаться от о2т, полученной при повторных наблюдениях 'для любого другого набора независимой переменной ут. Это требование также не всегда выполняется для реального эксперимента.
3. Независимые переменные x1, x2
x
I
!
1
it
0
и
is
1
to
I
та
i
и *
I
S
СО
0
1 1
ситуация 1. Получено уравнение регрессии первого порядка
Удалось настолько «хорошо» спланировать эксперимент, что получено уравнение регрессии вида (1) только первого порядка, но обладающее высокой степенью доверия и допустимой точностью. В этом случае показатели степени переменных хп равны 1, п = 1, 2,..., N. Оценка «хорошо» относится к тому обстоятельству, что был составлен так называемый ортогональный план эксперимента3, в результате чего исчезли слагаемые с сомножителями хпхи п, и = 1, 2, ..., N. Это означает, что получилось уравнение многомерной плоскости, на которой нет и не может быть экстремальных точек:
N
У = ь0 + ХЬпХп. (2)
п=1
Следует учесть, что в экономических задачах, как правило, существуют линейные ограничения вида
измеряются с пренебрежимо малой ошибкой по сравнению с ошибкой в определении координат У = f(X) или скаляра у.
При таких исходных предпосылках оказывается возможным вычислить коэффициенты Ь0 , ьп , Ьпи , Ьпп , ..., а также оценить их точность и точность уравнения регрессии (1) в целом.
Для того чтобы знать, что ждет экспериментатора по завершении имитационного эксперимента, сначала рассмотрим финишную часть экстремального эксперимента с имитационной моделью. Результатом могут быть три различные ситуации, которые различаются результатами регрессионного анализа:
— ситуация 1: получено уравнение регрессии первого порядка;
— ситуация 2: удалось получить уравнение регрессии второго порядка;
— ситуация 3: не удается получить удовлетворительное в смысле достоверности уравнение регрессии, как первого, так и второго порядков.
У anrxn <br, r = 1, 2, ..., R
^^ n,r n r' iii
n=1
n = 1, 2, ..., N
(3)
lXn ^ 0;
В таких случаях наиболее эффективный поиск оптимума (минимума или максимума) при числе переменных-факторов более двух — это метод линейного программирования. На рисунке 3 приведен соответствующий пример поверхности отклика — двухко-ординатная плоскость вида
у = Ь 0 + Ь1 х1 + Ь2 х2
с пятью линейными ограничениями, которые являются линиями, представляющими проекции многоугольника.
Из рисунка видно, что в этом случае решением являются координаты точки А в одном из углов поверхности отклика:
ут|П = А( хГ,х2™).
3 Построение такого плана имитационного эксперимента для уравнений регрессии первого порядка подробно рассмотрено в [7].
82
№ 3 (45) 2013
Издержки ▲
Неизвестная поверхность отклика
min min Оптимум (Х-|, Х2 )
Рис. 3. Оптимум на линейной поверхности отклика
Следует заметить, что при «плохом», т. е. неортогональном планировании, применить линейное программирование было бы проблематично. В этих случаях применяются методы нелинейной оптимизации, некоторые из которых будут рассмотрены ниже.
ситуация 2. Удалось получить уравнение регрессии второго порядка
Экспертно и с помощью серии опытов сделан вывод о наличии на поверхности отклика экстремума (их может быть и несколько). Эксперимент спланирован так, что получено уравнение регрессии второго порядка: показатели степени переменных xn равны 1 или 2, п = 1, 2,..., N. Опять составлен так называемый ортогональный план эксперимента4, в результате чего в уравнении регрес-
x2nu, n, u = 1, 2, ..., N. Уравнение регрессии
имеет вид:
У = b 0 +Х bnxn + Х bn
, x n
(4)
при наличии тех же ограничений (3).
Соответствующий пример показан на рис. 4.
В таких случаях применяется градиентный метод [8] или его разновидности. Градиентный метод. Поверхность ((X),
X = (x ■ x2
xN), на которой нужно наити
экстремальную точку, характеризуется градиентом — вектором, проекции которого являются частными производными по соответствующим переменным:
Vf (X) =
Э f( X)_ э f( X у _ э f( X)
Э x. Э x2
Э Xn
(5)
4 Построение такого плана для уравнений регрессии второго порядка также подробно рассмотрено в [7].
М "л2 "
В произвольной точке Хк, относящейся к области допустимых значений, градиент перпендикулярен касательной к изолинии
83
со
0
1 ¡1
сии нет слагаемых с сомножителями x x и
n=1
№ 3 (45) 2013
Рис. 4. Выпуклая вниз поверхность отклика
Область допустимых значений Хо.
уровня функции и указывает наиболее «крутое» направление изменения У(( X) из точ-
ки Xk.
Требуется: найти максимальное значение выпуклой вверх функции
f (Xi, x2
xs,
(6)
в общем случае — при отсутствии ограни чений:
< xn < те, n = 1, 2,
N.
Рис. 5. Градиентный метод поиска экстремума
Решение (классический вариант). В исходной точке из области допустимых значений Хк =^Х0 делается шаг по направлению
градиента в точку Хк, где значение функции явно меньше. Далее переходят к следующей итерации: к = к +1 и так далее, постепенно уменьшая шаг по мере приближения к экстремальной точке. Иллюстративно это можно представить в двухмерном виде (рис. 5).
84
№ 3 (45) 2013
Алгоритм прост, однако в известной степени требует от исследователя интуиции. Классический вариант решения используется в технике, в заводских лабораториях, но в экономике, как правило, не применяется.
Градиентный метод и его разновидности имеют два общих недостатка:
1) решения гарантированы, если все частные производные — непрерывные функции — как минимум на пути Х0, Х1, Х2, Х3, Х4 (рис. 5);
2) если экстремумов несколько, то для поиска всех экстремумов нужно искать новые «стартовые точки», т. е. для каждого экстремума должна быть определена своя точка Х0.
В экономических приложениях при использовании имитационного моделирования обычно имеют место ограничения вида (3), а иногда и дополнительные условия, например: некоторые или все х1 > 0, I = 1,2, ..., п могут принимать только целочисленные значения.
Разновидность градиентного метода, имеющая ограничения вида (3), называется методом Франка - Вульфа [1, 9]. В данном методе особенностью задачи поиска экстремума является то, что система ограничений содержит только линейные неравенства.
Эта особенность является основой для замены линейной функцией в окрестности исследуемой точки нелинейной целевой функции (чаще всего — неизвестной), благодаря чему решение исходной задачи сводится к последовательному решению задач линейного программирования. Казалось бы, это не упрощает решение задачи поиска экстремальной точки. Однако из приведенного ниже примера будет видно, что, если в результате эксперимента неизвестная поверхность отклика сведена к уравнению регрессии второго порядка, а размерность задачи невелика, то методы линейного программирования не потребуются.
Метод Франка - Вульфа. Процесс нахождения решения задачи начинается с оп-
ределения точки, принадлежащей области допустимых решений задачи. Пусть это точка Хк, тогда в ней вычисляется градиент функции (6) с помощью формулы (5):
Vf(Xk) =
д f(X,), Э f(X,), , Э f(X,)
д x1 ' д x2
д Xn
и строится линейная функция
F (Xk)=
д f( X,) Э f( Xk)
д X,
д x0
д f (X,
д X
X. (7)
Затем находится максимальное значение этой функции при ограничениях (3) путем определения особой точки Т.к.. Тогда за новое допустимое решение исходной задачи принимаются координаты точки Хк+1:
Хк+1 = хк + Хк(4 - хк), (8)
где Хк — некоторое число, называемое шагом вычислений, заключенное между нулем и единицей (0 < Хк < 1).
Шаг выбирается следующим образом:
К=
произвольное 0<Хк <1; обеспечиващее ma x{{( X^)} (9)
вточке X,
k+1
т. е. число Хк берется произвольно или определяется таким образом, чтобы значение функции f(Хк+1), зависящее от Хк, в точке Хк+1 было максимальным. Для этого необходимо найти решение уравнения
d f (k+Дк))
d Хк
= 0,
поскольку Хк+1 = Хк+1(Хк) — функция от Хк, и выбрать его наименьший корень. Если его значение больше единицы, то следует положить Хк = 1.
После определения числа Хк находят координаты точки Хк+1 и выясняется необходимость перехода к новой точке Хк+1 как к исходной для следующей итерации. Для этого исследуется необходимость проведения дальнейших вычислений с точностью 8, определяемой выражением
1 ¡1
85
к
5= f(X,+1) - f(X,) < е,
где е — допустимая погрешность.
(10)
I
!
1
it
0
и
is
1
to
I
та
i
*
I
S
t i
Если неравенство (10) не выполнено, то такая необходимость имеется, и производится переход к следующей итерации: к = к +1, опять выполняется этап 1 (см. приведенный ниже алгоритм), и т. д.
После конечного числа шагов и при выполнении неравенства (10) получается решение исходной задачи с необходимой точностью.
Итак, процедура нахождения экстремума функции (6) при ограничениях (3) градиентным методом состоит из следующих этапов:
Этап 1. Определятся исходная точка Хк = Х0 в области допустимых значений
для последующего решения задачи.
Этап 2. Находится градиент функции (6) в точке допустимого решения.
Этап 3. Строится функция (7) и находят ее максимальное значение при ограничениях (3).
Этап 4. Начинается вычисление новой точки Хк+1 по формуле (8).
Этап 5. Определяется шаг вычислений (9).
Этап 6. По формулам (8) находятся компоненты новой допустимой точки Хк+1.
Этап 7. Определяется точность 8 и проверяется выполнимость условия (10). Если условие не выполняется, то переходят к следующей итерации — опять к этапу 2, сделав замену к = к +1. Но если оно выполнено, то это означает, что найденная точка Хк+1 является приемлемым в смысле погрешности решением задачи.
Пример. Необходимо найти максимальное значение выпуклой квадратичной функции, которая в данном случае представляет собой уравнение регрессии второго порядка, полученного в результате эксперимента с имитационной моделью, выполненного по ортогональному плану.
f(X) = 2х1 + 4х2 - х2 - 2х22, (11) где X = (х1; х2) при условиях
х1 + 2 х2 < 8 2х1 - х 2 , х1 > 0, х 2 > 0
(12)
с точностью, определяемой неравенством5 Хк+1) - f(Хк)| < е, где е = 0,01 — допустимая погрешность.
Решение. Используем метод Франка -Вульфа. Найдем градиент функции
Vf (X) =
fd f(X) ; Э f (X)л д х1 ' д х2
= (2 - 2х1; 4 - 4х2
В качестве исходной точки из области допустимых значений возьмем Х0 = (0; 0).
Шаг 1. В точке Х0 = (0; 0) градиент целевой функции Vf( Х0) = (2; 4). Находим максимум функции
F (X0) = 2 х1 + 4 Х2
при ограничениях (12). Решением, обеспечивающим максимум6, является точка 4 = (0; 4).
Найдем следующую точку из области допустимых значений Х1 исходной задачи по формуле (8):
Х1 = Х0 + М4 - X0),
где 0 < ^ < 1.
Подставив вместо Х0 и 70 их значения, получим
у Хч = 0 • 0
1 п ^ „ . (13)
х2 = 0 + л0 • 4
Далее определим число Подставив в равенство (11) в качестве значений х1 и х2 их значения, полученные из соотношений (13), получим
5 Этот классический учебный пример нелинейного программирования без какой-либо привязки к имитационному моделированию приводился в различных учебниках, однако впервые его привел Иван Людвигович Акулич в [1] (Латвийский университет, Белорусский государственный экономический университет).
6 Решение в данном случае очевидно. Однако в более сложных случаях при наличии подобных линейных ограничений максимум придется искать, применяя симплекс-метод.
№ 3 (45) 2013
ДХ0) = -32Х0 + 16Х0,
найдем производную этой функции по Х0 и приравняем ее нулю:
б f(ХД0))
d Хо
= 16 - 64Х0 = 0,
откуда получим Х0 = 1/4. Поскольку 0 < Хк < 1, принимаем Х0 = 0,25 за величину шага 1.
В результате из формулы (13) получена точка Х1 = (0; 1). Осталось сделать проверку на точность: из формулы (11) получим f(X0) = 0, f(Х1) = 2, в результате - ЦХ0)\ = 2 > е = 0,01.
Вывод: решение не найдено.
Шаг 2. В точке Х1 = (0; 1) градиент целевой функции Vf( Х1) = (2; 0). Находим максимум функции F(X1) = 2х1 при ограничениях (12). Очевидным решением является точка
= (6,4; 0,8).
Найдем следующую точку из области допустимых значений Х2 исходной задачи по формуле (8): Х2 = Х1 + Х1( 71 - Х1). Это соотношение перепишем в виде
x1 = 6,4Х1 x2 = 1 - 0,8Х1
d f (Х2<Х1>)
= 12,8 - 83,152 Х1 = 0,
Шаг 3. В точке Х2 = (0,96; 0,97) градиент | целевой функции ЩХ2) = (0,08; 0,12). Нахо- 1 дим максимум функции F (Х2) = 0,08 х1 + 0,12 х2 ,2 при ограничениях (12). Решением является ^ точка г2 = (6; 0). ^
Найдем следующую точку из области допустимых значений Х3 исходной задачи по формуле (61): Х3 = Х2 + Х2(72 - Х2). Это соотношение перепишем в виде
x1 = 0,96 + Х2(6 - 0,96) = 0,96 + 5,04 Х2 x2 = 0,97 + Х2(0 - 0,97) = 0,97 - 0,97 Х2
. (15)
Далее, подставив в равенство (11) в качестве значений х1 и х2 их значения, полученные из соотношений (15), получим
f(Х2) = 2,9384 + 0,4032 Х2 - 27,3416 Х2
d f ((3 (Х2))
= 0,4032 - 54,6832 Х1 = 0.
(14)
Подставив в равенство (11) в качестве значений х1 и х2 их значения, полученные из соотношений (14), получим
f(Х1)=2+12,8Х1 -41,76Х2.
Далее найдем производную этой функции по Х1 и приравняем ее нулю:
б Х1
откуда следует, что Х1 = 0,153256... ~ 0,15.
Принимаем это значение за величину шага 2. Подставив в (14) значение Х1 = 0,15, получим: Х2 = (0,96; 0,97), f(X2) = 2,9966 и точность
Х2) - f( X)=2,9966- 2=0,9966 > е=0,01.
Вывод: опять решение не найдено, но мы приблизились к нему.
б Х2
Определяем шаг: Х2 = 0,007373... « 0,007, и в результате получим третью точку
Х3 = (0,99528; 0,96321), f(X3)=2,99957=2,9996, и точность
^ (Х3) - f (Х2)| = 2,9996 - 0,9966 = = 0,003 < е = 0,01.
Р е ш е н и е н а й д е н о : т о ч к а Х3 = (0,99528; 0,96321) является искомым решением поставленной задачи. Задав меньшую величину е, можно было бы, после дополнительных шагов, еще ближе подойти к точке максимального значения целевой функции.
ситуация 3. не удается получить удовлетворительное уравнение регрессии, как первого, так и второго порядков
Обстоятельства таковы, что экспериментатору по каким-то причинам не удалось составить ортогональный план (например, он не владеет соответствующей методикой). Ситуация нередкая, даже, можно сказать, — типичная. Тем не менее, и в ней есть
87
I
!
Ü
Ü о
и
is
К
to
I
та
i
I
I
s
to t l
возможность оперативного планирования и управления экспериментом, используя разновидность градиентного метода — метод покоординатного спуска Гаусса - Зей-деля [2, 7]. Этот метод при проведении имитационных экспериментов проще в реализации, чем классический градиентный метод, за счет того, что на очередной итерации спуск осуществляется постепенно вдоль каждой из координат. Кроме того, в условиях, когда вид функции неизвестен, об аналитическом вычислении градиента можно забыть. Вся информация может быть получена только из опыта. Основное достоинство метода заключается в том, что как вид функции f(x1, х2, ..., хм), так и ее градиент Vf( Хк) могут оставаться неизвестными во всех точках Хк, где к — номер шага в процедуре поиска экстремума.
Метод покоординатного спуска Гаусса - Зейделя. Сущность алгоритма покоординатного спуска очень простая. Допустим, что нужно найти минимум на сложной поверхности f( х1, х2) на рис. 6. Введем обозначения для точек Х0, Х1, Х2, Х3, Х4:
X0 = а; Х1 = Ь; Х2 = с; Х3 = d; Х4 = е.
Шаг 1. Выбираем произвольную точку Х0 =а из области допустимых значений
Область допустимых значений
Рис. 6. Пошаговый спуск по методу Гаусса - Зейделя
с координатами {x(a), xif}. Фиксируем значение {x<a)} и ищем минимум функции
min f (x1, x«a)).
xi
Поскольку аналитическая формула f (x1, x2) неизвестна, то поиск точки делается так: простейшим известным «методом деления пополам», описанной в учебной литературе по численным методам, с помощью серии опытов точка X2 = b определяется с заданной допустимой погрешностью А — величиной, с которой сравнивается последний использованный отрезок при реализации метода деления пополам. Причем, если есть дополнительное ограничение для x1, в соответствии с которым этот фактор может принимать только целочисленные значения, то количество опытов для определения X2 = b сокращается.
Очевидно, что x2a) = x2b).
Шаг 2. В найденной в области допустимых значений точке X1 = b с координатами {xjb), x2b)} фиксируем значение {x(b)} и ищем минимум функции
min f (x(b), x2).
x2
Поиск точки опять делается «методом деления пополам», причем точка X2 = c определяется также с допустимой заданной погрешностью.
Заметим, что x(b) = xjc).
Шаг 3. В найденной точке X2 = c с координатами {x(c), x2c)} фиксируем значение {x2c)} и ищем минимум функции
min f (x1, x2c))
xi
«методом деления пополам». Определяем положение точки X3 = d с заданной погрешностью.
Очевидно, что x2c) = x2d).
Шаг 4. В X3 = d с координатами {x(d), xif)} фиксируем значение {x(d)} и ищем минимум функции
Ad)
№ 3 (45) 2013
Поиск точки Х4 = е опять делается методом деления пополам с допустимой погрешностью.
Отметим, что х(б) = х\е).
Критерий останова. После каждого шага исследуется необходимость проведения дальнейших вычислений с точностью 8, определяемой неравенством
5 = ^ (Хк+1) - f (Хк)| < е,
где е — допустимая погрешность определения экстремального значения функции f(xv х2); Хк+1 — последняя точка; Хк — предпоследняя. Если это неравенство не выполнено, то такая необходимость имеется, и производится переход к следующему шагу, и т. д.
Замечание. Этот метод хорошо работает и в следующих случаях:
— известна аналитическая формула f (х„ хг);
— с помощью ортогонального планирования удается получить удовлетворительное в смысле достоверности уравнения регрессии второго порядка.
В этих случаях также не нужно определять сам градиент. Рассмотрим пошаговое решение, которое годится для обоих этих случаев.
Шаг 1. В точке Х0 = а фиксируем значение {х2а)}. Ищем минимум функции
min f х«а)).
Для этого найдем первую и вторую частные производные f (х1, х2). Делаем шаг из точки Х0 = а в точку Х2 = Ь, где частные производные удовлетворяют соотношениям
—f( Э х1 1
x2a)
) = 0
/ V -Ub)
Э2
d (xi)2
-f(х1,x2a)) > 0.
V 1 2 / v_v(ь)
Отметим, что x2a) _ x2b).
Шаг 2. В точке X1 = b с координатами [x'f, x«b)} фиксируем значение {xf} и ищем минимум функции
min f (x\b\ x2).
Опять классическим способом находится точка Х2 = с с помощью первой и второй частных производных, причем в найденной точке Х2 = Ь частные производные
д
Э x Э2
f (x(b),x2 )x2 _x,0
f ( x(b) ) >0;
д (Х2)2
Заметим, что х2а) = х2Ь).
И так далее. Критерий останова тот же самый.
Метод Гаусса - Зейделя обладает теми же недостатками, что и градиентный метод, особенно в случае «острого» оврага, где производные имеют разрыв.
Теорема о сходимости метода покоординатного спуска. Необходима особая интуиция для выбора начальной точки, поскольку из одной начальной точки можно попасть только в одну экстремальную точку. А как быть в случае, когда их несколько?
Для простоты рассмотрим функцию двух переменных f(x1, х2). Выберем некую стартовую точку (х(0), х20)) и проведем изолинию уровня через эту точку. Выбор начальной точки упрощается благодаря следующей теореме.
Теорема (без доказательства). Пусть в области G, ограниченной этой линией уровня, выполняются неравенства, означающие положительную определенность квадратичной формы:
д2 д2
f(xvх2) > а > 0 , ^f(xvх2) > Ь > 0, Х1 Х2
--^(х1,х2)| < с , а• Ь > с2.
Х1 Х2
Тогда спуск по координатам сходится к минимуму из начального приближения
(x!
(0) v(0) 1
xf), причем линеино.
Благодаря простоте программной реализации метод Гаусса - Зейделя пользуется популярностью, особенно у молодых исследователей, судя по публикациям в российской периодике за последнее десятилетие.
I
Uj
89
и
№ 3 (45) 2013
Заключение
В настоящей статье автор обратил внимание на то, что к опыту нужно тщательно подготовить имитационную модель, проверить исходные данные и продумать, как в модели будут обрабатываться результаты. Тщательная подготовка будет способствовать уменьшению ошибки опыта, которая является суммарной величиной, состоящей из ряда ошибок: ошибок при измерении факторов, параметра оптимизации и ошибок при проведении опыта. Ошибки подразделяются на случайные и систематические. Для того чтобы компенсировать влияние систематических ошибок, опыты нужно рандомизировать во времени. Тщательное, скрупулезное выполнение эксперимента, несомненно, является главным условием успеха исследования.
Рассмотренные выше методы позволяют упорядочивать серии опытов, которые иногда напоминают «случайные блуждания» в области допустимых значений, и сделать эксперимент конечным по времени, причем в результате планирования время эксперимента ^ для решения задачи поддается расчету. | В литературе по применению акторного Ц имитационного моделирования приводится | много практических примеров с решениями, | однако ответы на вопросы о том, как были | получены рациональные решения, как пра-§ вило, отсутствуют, так как основное внима-
I ние уделялось теории и практике имитаци-* онного моделирования. Чтобы снять этот § вопрос, следует отметить, что, например, ^ методом Гаусса - Зейделя выполнена автоматически оптимизация решений в таких
| РНдпт-моделях, как [7, глава 10]:
« — замкнутая модель «Минимизация за-
Ц трат на швейной фабрике»;
| — двухслойная модель «Финансовая ди-
II намика компании» для адаптивного выбора [з вариантов бизнес-процесса.
Однако в более сложных случаях на вто-| ром месте по распространенности исследо-§ ватели проводят ортогональное планирова-| ние, получая уравнения регрессии второго § порядка в области экстремального значе-
ния неизвестной поверхности отклика, с последующим применением метода Франка-Вульфа.
Уравнения регрессии первого порядка используются в простейших случаях.
Список литературы
1. Акулич И. Л. Математическое программирование в примерах и задачах. М.: Высшая школа, 1986. — 320 с.
2. Бодров В. И., Лазарева Т. Я., Мартемья-нов Ю. Ф. Математические методы принятия решений. Тамбов: ТГТУ, 2004. — 124 с.
3. Босс В. Лекции по математике. Т. 4. Вероятность, информация, статистика. М.: КомКнига, 2005. — 216 с.
4. Девятков В. В. Мир имитационного моделирования: взгляд из России // Прикладная информатика. 2011. № 4 (34). С. 9-29.
5. Емельянов А. А. Концепция и возможности ак-торно-ориентированной системы имитационного моделирования Actor Pilgrim // Прикладная информатика. Ч. I — 2012. № 6 (42). С. 49-66; Ч. II — 2013. № 1 (43). С. 41-53.
6. Емельянов А. А. Лаг-генераторы для моделирования рисковых ситуаций в системе Actor Pilgrim // Прикладная информатика. 2011. № 5 (35). С. 98-117.
7. Емельянов А. А., Власова Е. А., Дума Р. В., Емельянова Н. З. Компьютерная имитация экономических процессов / под ред. А. А. Емельянова. М.: Маркет ДС, 2010. — 464 с.
8. Коршунов Ю. И. Математические основы кибернетики. М.: Энергия, 1980. — 424 с.
9. Мицель А. А, Шелестов А. А. Методы оптимизации. Ч. 1. Томск: ТУСУР, 2002. — 192 с.
10. Монтгомери Д. К. Планирование эксперимента и анализ данных. Л.: Судостроение, 1980. — 384 с.
11. Налимов В. В., Чернова Н. А. Статистические методы планирования экстремальных экспериментов. М.: Наука, 1965. — 342 с.
12. Тейлор Дж. Введение в теорию ошибок. М.: Мир, 1985. — 272 с.
13. Терентьев Н. Е. Риски инновационного развития и повышение конкурентоспособности компании // Современная конкуренция. 2008. № 3 (9). С. 69-77.
