Известия РАН. Энергетика, 2021, № 4, стр. 24-41

ВВЕДЕНИЕ

Основная задача анализа балансовой надежности (БН) электроэнергетической системы (ЭЭС) – оценить вероятностные показатели (вероятность, частота, математическое ожидание (МО) и др.) дефицита мощности (ДМ) в ЭЭС, наличие которого рассматривается как нарушение (отказ) нормального функционирования ЭЭС [1–3]. Аналогом упомянутых показателей балансовой надежности за рубежом являются: LOLP (Loss of load probability), LOLF (Loss of load frequency), EPNS (Expected power not supplied) [4]. Их сравнительная оценка, соответствие и формулы преобразования описаны в [5].

Случайное состояние системы, в основном, характеризуется вероятностным характером нагрузки $L = \left\{ {{{L}_{i}},~i = 1, \ldots ,n} \right\},$ располагаемой генерации $G = \left\{ {{{G}_{i}},~i = 1, \ldots ,n} \right\}$ и топологии электрической сети. При этом дефициты мощности могут иметь как локальный (в отдельных подсистемах), так и глобальный (по ЭЭС в целом) характер. Локальные ДМ определяются, в основном, состоянием (включен – отключен) и пропускной способностью элементов электрической сети.

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

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

На втором этапе ЭЭС рассматривается как система с ограниченной ПС МСС (структура 2). Здесь МСС с ограниченной ПС объединяют некоторое множество концентрированных ЭЭС, которые в России принято определять как зоны надежности [7] или зоны свободного перетока мощности [8]. Основной направленностью задач БН здесь является анализ возможности появления локальных дефицитов мощности, вызванных ограниченностью ПС МСС. Именно такая структура является объектом анализа в представленной работе.

Полная детализация ЭЭС с учетом систем управления режимами является предметом анализа режимной надежности; ориентирована на задачи развития систем автоматики и диспетчерского управления и в данной работе не рассматривается.

Следует отметить, что все технические системы, в том числе и ЭЭС, характеризуются малой (доли процента) вероятностью событий типа “отказ” (в БН ЭЭС – это дефицитное состояние ЭЭС). При этом, как правило, отказ проявляется при бесконечно большой комбинации внешних событий. Так в ЭЭС с двумя узлами нагрузки и фиксированной предельной мощностью электрических станций дефицит мощности может проявляться при недопустимо большом повышении нагрузки либо первого, либо второго, либо того и другого узлов. При этом величина превышения суммарной нагрузки над генерацией является непрерывной и отсюда бесконечно неопределенной в реализации величиной. Число узлов нагрузки и генерации реальных расчетных схем ЭЭС измеряется тысячами. Отсюда многократно возрастает неопределенность дефицитного состояния ЭЭС, становится затруднительным процесс идентификации наиболее значимых комбинаций событий, приводящих к отказу системы.

Основным математическим аппаратом при анализе надежности ЭЭС в настоящее время являются методы статистических испытаний и, в частности, метод Монте-Карло (ММК) [9–13 и др.]. Данные методы позволяют моделировать системы со сложными функциональными взаимосвязями, не поддающимися аналитическому описанию, в том числе учитывать стохастическую неопределенность, связанную с предложением, спросом и пропускной способностью межсистемных связей.

Известен основной недостаток ММК – для обеспечения точности результатов требуется достаточно большое число $N$ испытаний, поскольку $N$ обратно пропорционально вероятности моделируемых событий. В частности, при моделировании события с вероятностью 10^–3 и погрешностью моделирования 1% требуется не менее 10⁵ статистических испытаний [14]. При анализе балансовой надежности реальных ЭЭС наиболее значимыми являются кратные отказы (одновременный отказ двух линий электропередачи (ЛЭП), одновременный отказ ЛЭП и энергоблока, отказ энергоблока в период пиковых нагрузок и др.), поскольку единичные отказы электрооборудования ЭЭС как правило не приводят к ограничению электропотребления – ЭЭС проектируются и сооружаются с учетом критерия N-1, при котором единичный отказ любого элемента ЭЭС не должен приводить к отказу процесса электроснабжения потребителей. Вероятность кратных независимых отказов равна произведению вероятностей отказа отказавших элементов. В результате вероятность кратных отказов принимает значения 10^–6–10^–8. При этом для идентификации, моделирования и анализа послеаварийных состояний ЭЭС методом ММК требуется число испытаний не менее 10⁷–10⁹. При меньшем числе испытаний идентификация таких редких событий становится маловероятной, что приводит к значимой погрешности результирующих показателей БН. С учетом того, что при каждом испытании требуется выполнение оптимизационных расчетов (оптимальное потокораспределение) анализ БН ЭЭС потребует несколько часов машинного времени. Следует отметить, что расчет показателей надежности, как правило, является составной частью более общих задач, например, выбор оптимальной конфигурации электрической сети при ее развитии. Отсюда длительность расчета надежности варианта развития ЭЭС как отдельного расчетного блока лимитируется по времени. Современная энергосистема требует все более сложной детализации и расширения учитываемых факторов, что нелинейно увеличивает сложность задачи оценки надежности ЭЭС. В совокупности это приводит к нецелесообразности, а иногда и невозможности использования стандартного ММК для оценки БН энергосистемы. Требуются модификации ММК, направленные на увеличение вычислительной эффективности статистических методов.

Одним из путей решения проблемы выполнения большого числа относительно однотипных расчетов, связанных с моделированием в многомерном вероятностном пространстве редких (с вероятностью отказов менее 10^–5) событий является идея параллельных вычислений. Технология одновременного использования нескольких компьютеров (или процессоров) активно развивается в течение последних 30 лет [15–17]. Параллельные вычисления, безусловно, снижают вычислительную нагрузку ММК. Однако зависимость между количеством параллельных вычислительных ресурсов и временем вычислений близка к линейной функции [17]. Отсюда положительный эффект распараллеливания вычислений становится не столь очевидным при моделировании редких событий.

Другим направлением является использование методов паттернов, метамоделей, искусственного иммунитета, искусственной нейронной сети, опорных векторов наименьших квадратов [18–21] и др. Эти методы также показали снижение вычислительной нагрузки метода ММК. Однако, как правило, для задач БН ЭЭС они не апробированы, ориентированы на определенный класс задач, зачастую малоэффективны и не обеспечивают требуемую точность и надежность идентификации событий с вероятностью менее 10^–4, что характерно для реальных систем энергетики. Отсюда проблема идентификации редких событий остается по-прежнему актуальной.

Более универсальный характер имеют итерационные методы, связанные с преобразованием анализируемых пространств, функций вероятностных распределений случайных переменных и критериальных функций. Наиболее распространенными в данном классе методов являются методы подпространств и значимой выборки [22–25], а также кросс-энтропийные методы [26–28].

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

ДЕФИЦИТ МОЩНОСТИ КАК РЕДКОЕ СОБЫТИЕ

Пусть $x \in {{\mathbb{R}}^{m}}$ – случайный вектор, который объединяет все вероятностные входные переменные ${{L}_{i}},{{G}_{i}},$ и др. В простейшем случае, когда учитываются только располагаемая генерация и нагрузка и не учитываются случайные состояния элементов системы, температуры окружающей среды и др., $x = \left\{ {{{L}_{i}},{{G}_{i}},~~i = 1, \ldots ,n} \right\},$ где $i~$ – номер узла электрической сети. С целью концентрации внимания на механизмах идентификации редких событий, связанных с отказами функционирования ЭЭС множество случайных, вероятностно определенных величин в данной работе ограничено только нагрузкой и располагаемой генерацией.

Дефицитность системы определяется разностью $~{{D}_{{{\Sigma }}}} = ~{{L}_{{{\Sigma }}}} - {{G}_{{{\Sigma }}}}$ = $\sum {{{L}_{i}}} - \sum {{{G}_{i}}} ,\,\,{{L}_{{{\Sigma }}}} > {{G}_{{{\Sigma }}}}.$ Ее можно оценить, например, с$~$ помощью функции избыточности системы, $S = {{G}_{{{\Sigma }}}} - {{L}_{{{\Sigma }}}}$ или (при учете пропускной способности электрической сети) $S = \min \left( {{{G}_{i}} - {{L}_{i}},~\forall i} \right).$ При независимости $\left\{ {{{L}_{i}},{{G}_{i}}} \right\}$ часто имеет смысл рассматривать обобщенную случайную величину – располагаемая генерирующая мощность узла, ${{r}_{i}} = {{G}_{i}} - {{L}_{i}}$. Ее МО и дисперсия: ${{\mu }_{{{{r}_{i}}}}} = {{\mu }_{{{{G}_{i}}}}} - {{\mu }_{{{{L}_{i}}}}};$ $\sigma _{{{{r}_{i}}}}^{2} = \sigma _{{{{G}_{i}}}}^{2} + \sigma _{{{{L}_{i}}}}^{2}.$ В результате число управляющих стохастических переменных сокращается до числа узлов электрической сети. При этом ${{D}_{{{\Sigma }}}} = - {{r}_{{{\Sigma }}}} = - \sum {{{r}_{i}}} ,\,\,\,\,{{r}_{{{\Sigma }}}} < 0.$

В общем случае при анализе локальных ДМ необходимо учитывать законы распределения мощности в электрической сети. Каждый узел характеризуется экспортом ${{u}_{i}},~i = 1, \ldots ,n$ или импортом (–${{u}_{i}}$) мощности, определяемыми диспетчерским управлением ЭЭС и зависящими от располагаемой мощности узлов ${{u}_{i}} = {{u}_{i}}\left( r \right).$ При этом локальный ДМ ${{D}_{i}} = {{u}_{i}} - {{r}_{i}} = {{u}_{i}} - {{G}_{i}} + {{L}_{i}};$ ${{D}_{i}} > 0.$ При локальных ДМ системный ДМ фиксируется, если хотя бы в одном из узлов имеет место локальный ДМ $~{{D}_{{{\Sigma }}}} = \sum {{{D}_{i}}} .$ Для более точного учета дефицитности ЭЭС необходимо учитывать потери активной мощности в сети, которые определяются распределением потоков мощности $z = \left\{ {{{z}_{j}},j = 1, \ldots ,m} \right\}$ в элементах электрической сети и имеют нелинейный (квадратичный) характер. Учет потерь мощности связан с появлением дополнительного условия ${{\pi }_{{{\Sigma }}}} = \sum {{{u}_{i}}} .$ При этом ${{D}_{{{\Sigma }}}} = ~{{L}_{{{\Sigma }}}} - {{G}_{{{\Sigma }}}} + ~{{\pi }_{{{\Sigma }}}},$ ${{D}_{{{\Sigma }}}} \geqslant 0.$

Вектор $z$ определяется решением задачи оптимального распределения мощности нагрузки между источниками питания. В простейшем случае, приемлемом в задачах оценки надежности, можно считать $z = Cu{\text{,}}$ где $C$ – матрица потокораспределения [7]. В результате ДМ в ЭЭС представляет сложную (часто определяемую алгоритмически) функциональную зависимость от совокупности управляющих переменных, имеющих случайный характер.

Состояние дефицитности системы фиксируется, если некоторая заданная на множестве управляющих переменных критериальная функция $\psi (x)$ становится меньше заданного предельного значения ${{\psi }_{{{\text{lim}}}}}{\text{:}}$ $\psi \left( x \right) < {{\psi }_{{{\text{lim}}}}}.$ В практике для идентификации дефицитности чаще используется функция с нулевым порогом, $S\left( x \right) = \psi \left( x \right) - {{\psi }_{{{\text{lim}}}}},$ которая при дефицитности системы принимает, например, отрицательное значение, $S\left( x \right) < 0.$ При этом область

определяет совокупность дефицитных состояний системы.

При учете локальных ДМ в ЭЭС критериальная функция имеет вид:

Критериальной функции соответствует индикаторная функция$:$

которая позволяет выразить многие логические конструкции в аналитическом виде. В частности, вероятность ${{\mathbb{P}}_{H}}$ состояния отказа (наличие ДМ в ЭЭС) может быть определена как МО индикаторной функции $J\left( x \right),$ принимающей в состоянии ЭЭС $x$ значение 1 при наличии отказа в ЭЭС и 0 при его отсутствии: где $f\left( x \right)$ – плотность распределения (ПР) определяющей рассматриваемое событие многомерной случайной величины $x.$

При заданной ПР $f\left( x \right)$ вероятность появления редкого события определяется областью H, удовлетворяющую критерию $S\left( x \right) < 0.$

Функцию $F\left( x \right)$ (или плотность$~f\left( x \right)$) распределения случайного вектора $x$ состояний системы практически невозможно представить в виде аналитического выражения, поскольку она является результатом свертки разнотипных, с различающимися параметрами большого числа параметров системы (нагрузки, генерации, бинарные переменные, определяющие состояние элементов системы, солнечная радиация, скорость ветра и др.). Именно поэтому основным методом анализа надежности ЭЭС в настоящее время является метод статистических испытаний, позволяющий путем многократных испытаний и их статистической обработки получить достаточно точную оценку искомых величин. При этом в качестве исходных данных, как правило, задаются параметры маржинальных распределений исходных случайных величин (например, нагрузка узла $i$ ЭЭС описывается нормальным (Гауссовским) распределением с математическим ожиданием ${{\mu }_{i}}$ и среднеквадратическим отклонением (СКО) ${{\sigma }_{i}},$ ${{L}_{i}}\sim N(x,{{\mu }_{i}},~{{\sigma }_{i}}).$

МЕТОД ПОДПРОСТРАНСТВ СОСТОЯНИЙ

Эффективность метода Монте-Карло может быть повышена за счет применения методов уменьшения дисперсии, таких как выборка по значимости. Основная идея данной группы методов – формирование последовательности подмножеств ${{H}_{1}} \supset {{H}_{2}} \supset \ldots \supset {{H}_{m}} = H$ пространства состояний системы, где каждое последующее подпространство увеличивает вероятность идентификации редких событий и определяется на базе предыдущего, образуя последовательность цепи Маркова. При этом вероятность появления редкого события

Каждое последующее подмножество выбирается так, чтобы вероятность условного события $\Pr \left( {{{H}_{j}}{\text{|}}{{H}_{{j - 1}}}} \right)$ была бы достаточно велика. В результате малая вероятность представляется произведением относительно больших вероятностей [23–25].

Одним из путей формирования $\left( {{{H}_{j}}{\text{|}}{{H}_{{j - 1}}}} \right)$ является отбор заданной доли ${{p}_{0}}$ наиболее значимых событий $X_{b}^{j} = \left( {x_{1}^{j},~x_{2}^{j}, \ldots ,x_{k}^{j}} \right),$ где $k = {{p}_{0}}N;$ $N$ – объем выборки. Индекс b характеризует максимальный для множества $X_{b}^{{j - 1}}$ уровень критериальной функции $\psi \left( x \right).$ Его более детальное определение дано ниже). Значимость событий определяется по величине $\psi \left( x \right)$ – чем меньше $\psi \left( x \right),$ тем больше значимость $x$ (при $\psi \left( {{{x}_{i}},i = 1,...k} \right) < 0$ вся область $X_{b}^{j}$ состоит из значимых событий типа “отказ” (дефицитные состояния ЭЭС)). Дифференциация по значимости определяет механизм формирования множества $X_{b}^{j}.$ Полная, полученная на основе данных этапа $j - 1$ выборка ${{H}_{{j - 1}}}$ упорядочивается по возрастанию $\psi \left( x \right){\text{:}}$ $\psi \left( {{{x}_{1}}} \right) \leqslant \ldots \leqslant \psi \left( {{{x}_{N}}} \right).$ Первые ${{p}_{0}}N~$ событий (${{p}_{0}} \times 100$ перцентиль функции $\psi \left( x \right)$) определяют множество $X_{b}^{j}$ и соответствующую ему максимальную величину критериальной функции, ${{b}_{j}} = {\text{max}}\left( {\psi \left( x \right),~x \in X_{b}^{j}} \right),$ которая, в свою очередь, является основой для формирования нового множества ${{H}_{j}} = \left\{ {~x{\kern 1pt} :\,\,~\psi \left( x \right) < {{b}_{j}}} \right\}.$ При этом величину ${{p}_{0}}$ можно рассматривать как вероятность условного события $\Pr \left( {{{H}_{j}}{\text{|}}{{H}_{{j - 1}}}} \right) = {{p}_{0}}.$

Представленный поэтапный процесс формирования множества редких событий характеризуется положительной величиной ${{b}_{j}} > 0$ на всех промежуточных этапах. Это означает, что множество $X_{b}^{j}$ содержит как события-отказы, $\psi \left( x \right) < 0,$ так и не являющиеся отказами $~\psi \left( x \right) > 0,$ т.е. принцип отбора событий сводится к исключению менее значимых событий и расширению области более значимых событий, ${{b}_{j}} < {{b}_{{j - 1}}}.$ Изначально ${{b}_{0}} = \infty ,$ что означает, принадлежность зоне анализа всех сгенерированных на базе маржинальных ПР состояний системы. Однако на последующих этапах при генерации анализируемого множества псевдослучайных состояний системы вводится ограничение: $\psi \left( x \right) < {{b}_{{j - 1}}}.$ Как правило, это реализуется заменой параметров ПР случайных переменных некоторыми новыми расчетными величинами.

На последнем этапе ${{b}_{j}} < 0.$ Это означает, что все события в выделяемой по принципу усечения множества с долей ${{p}_{0}}$ являются отказами. Но отказами могут быть и события, не попавшие в определяемую перцентилем область, ${{b}_{j}} < \psi \left( x \right) < 0.$ Здесь условная вероятность определяется согласно соотношению: $\Pr \left( {{{H}_{m}}{\text{|}}{{H}_{{m - 1}}}} \right) = {{N\left( {\psi \left( x \right) < 0} \right)} \mathord{\left/ {\vphantom {{N\left( {\psi \left( x \right) < 0} \right)} N}} \right. \kern-0em} N},$ где $N\left( {\psi \left( x \right) < 0} \right)~$ – число элементов выборки объема N, удовлетворяющих требованию $\psi \left( x \right) < 0.$

В зависимости от алгоритма полный объем выборки ${{N}^{{(j)}}}$ на промежуточных этапах может отличаться от заданного N, однако при этом сохраняется принцип отбора – к дальнейшему рассмотрению принимается только ${{p}_{0}}$ – часть сформированного для анализа множества состояний системы.

АЛГОРИТМ АДАПТИВНОЙ ВЫБОРКИ

1. Инициализация: $k = 1;$ ${{N}_{c}} = {{p}_{0}}N;$ ${{N}_{s}} = {1 \mathord{\left/ {\vphantom {1 {{{p}_{0}}}}} \right. \kern-0em} {{{p}_{0}}}};$ $\lambda = 0.6.$

2. Генерация по стандартному нормальному (Гауссовскому) распределению $N$ псевдослучайных чисел: $U = \left\{ {{{u}_{1}}, \ldots ,{{u}_{N}}} \right\}.$

3. Преобразование множества $U$ в матрицу $X = \left\{ {{{x}_{1}}, \ldots ,{{x}_{N}}} \right\}$ именованных случайных переменных (нагрузка, генерация, состояние элементов системы), согласно их маржинальным законам распределения ${{x}_{i}} = {{{{\varphi }}}_{i}}\left( {{{u}_{i}}} \right).$

4. Определение вектора критериальных функций $\psi = \left\{ {{{\psi }_{i}}\left( {{{x}_{i}}} \right)} \right\}$ и его сортировка по возрастанию функции: ${{\psi }_{i}}\left( {{{x}_{i}}} \right) \geqslant {{\psi }_{{i - 1}}}\left( {{{x}_{{i - 1}}}} \right).$

5. Определение на множестве $\psi $ перцентиля $\left( {{{p}_{0}} \times 100} \right),$ соответствующему ему параметра ${{b}_{k}} \geqslant 0$ и множества именованных значимых состояний системы ${{H}_{{xk}}} = \left\{ {{{x}_{i}}:{{\psi }_{i}}\left( {{{x}_{i}}} \right) < {{b}_{k}}} \right\}.$

6. Генерация новой выборки N случайных состояний системы на базе полученных в п. 5 опорных состояний системы согласно алгоритму адаптивной выборки, где после определенного числа генераций меняются МО и СКО нормального распределения, согласно которому осуществляется генерация случайных чисел.

7. Пп. 4–6 повторяются до тех пор, пока ${{b}_{k}} \geqslant 0.$ При этом условная вероятность ${\text{Pr}}\left( {{{H}_{k}}\left| {{{H}_{{k - 1}}}} \right.} \right) = {{p}_{0}}.$ На последнем шаге m, при ${{b}_{m}} < 0,$ определяется число ${{N}_{m}}{\text{\;}}$элементов, удовлетворяющих условию ${{\psi }_{{{{N}_{m}}}}}\left( {{{x}_{{{{N}_{m}}}}}} \right) \leqslant 0,$ ${{\psi }_{{{{N}_{m}} + 1}}}\left( {{{x}_{{{{N}_{m}} + 1}}}} \right) > 0.$ Результирующая вероятность редкого события $\Pr \left( H \right) = {{p_{0}^{{m - 1}}{{N}_{m}}} \mathord{\left/ {\vphantom {{p_{0}^{{m - 1}}{{N}_{m}}} N}} \right. \kern-0em} N}.$

8. Определение остальных анализируемых вероятностных показателей (МО дефицита мощности и др.).

КРОСС-ЭНТРОПИЙНЫЙ МЕТОД

Кросс-энтропийный метод (КЭМ) [26–28] основан на замене реальной функции $f\left( x \right)$ плотности распределения анализируемой многомерной случайной величины, по которой делается выборка ММК, некоторой вспомогательной функцией плотности распределения (ВПР) $q\left( x \right),$ смещающей сферу анализа в область представляющего интерес редкого события. Данный подход широко известен при вычислении интегралов сложных функций статистическими методами, где подынтегральная функция умножается и делится на некоторую функцию ПР, полностью определенную на рассматриваемом интервале интегрирования:

где $W\left( x \right) = {{f\left( x \right)} \mathord{\left/ {\vphantom {{f\left( x \right)} {q\left( x \right)}}} \right. \kern-0em} {q\left( x \right)}};$ $\mathbb{E}\left[ \ldots \right]$ – оператор математического ожидания.

При этом возникает проблема выбора вспомогательной функции – необходимо, чтобы основанная на ней расчетная процедура была бы не только адекватной по направленности (смещения анализируемой области в сторону редкого события), но и эффективной по быстродействию и сходимости. При наличии ВПР вероятность дефицита мощности (2), может быть представлена в виде математического ожидания (МО) взвешенной индикаторной функции $J\left( x \right)W\left( x \right),$ определенной в многомерном пространстве случайных величин с плотностью распределения $q\left( x \right){\text{:}}$

где $W(x) = {{f(x)} \mathord{\left/ {\vphantom {{f(x)} {q(x)}}} \right. \kern-0em} {q(x)}}$ – весовая функция, определяемая исходной $(f(x))$ и вспомогательной $\left( {q\left( x \right)} \right)$ плотностями распределения.

Нетрудно видеть, что в выражении (3) $q\left( x \right)$ может быть любой функцией, лишь бы она удовлетворяла требования, предъявляемым к ПР. При этом не обязательно, чтобы она принадлежала к тому же классу ПР, что и $f\left( x \right).$ Как правило, в качестве $q\left( x \right)$ принимается ПР многомерного нормального распределения.

При моделировании случайных состояний системы $\{ {{x}_{i}},i = 1,....,N\} $ методом Монте-Карло на базе плотности распределения $q\left( x \right)$ и набора дискретных состояний системы $\left\{ {{{x}_{i}}} \right\}$ вероятность (3) редкого события трансформируется к виду:

где $\left\{ {W\left( {{{x}_{i}}} \right)} \right\}$ может рассматриваться как совокупность весовых коэффициентов.

Дисперсия $\mathbb{V}({{\hat {\mathbb{P}}}_{F}})$ оценки вероятности определяется как средняя дисперсий случайной величины, распределенной по закону Бернулли:

Коэффициент вариации $\mathbb{C}\mathbb{V}{\text{,}}$ как мера ошибки оценки измеряемой случайной величины (вероятности ${{\hat {\mathbb{P}}}_{H}}$), определяется отношением стандартного отклонения к МО. В частности, для квадрата $\mathbb{C}{{\mathbb{V}}^{2}}\left( {{{{\hat {\mathbb{P}}}}_{H}}} \right){\text{:}}$

Отсюда видно, что коэффициент вариации (КВ) при моделировании по методу Монте-Карло обратно пропорционален $\sqrt {N{{\mathbb{P}}_{F}}} $ и поэтому для получения достаточно точной оценки вероятности рассматриваемого события при малых ${{\mathbb{P}}_{F}}$ для обеспечения приемлемого $\mathbb{C}\mathbb{V} = 0.1 - 0.2$ необходимо чрезвычайно большое количество N испытаний ($\sim {{10}^{5}}$ при $\mathbb{C}\mathbb{V} = 0.1$ и ${{\mathbb{P}}_{F}} = 0.001$).

В качестве критерия оптимизации при выборе наиболее эффективной в классе $\left\{ {q\left( x \right)} \right\}$ функции $q{\text{*}}(x)$ может служить минимум дисперсии ${{\mathbb{P}}_{F}}{\text{:}}$

Теоретически наилучшей функцией, приводящей к нулевой дисперсии искомой оценки вероятности ${{\mathbb{P}}_{H}}$, является функция [28]:

Поскольку оптимальная ВПР зависит от неизвестных величин ${{\mathbb{P}}_{H}}{\text{\;}}$и $J\left( x \right)$ непосредственное аналитическое определение $q{\text{*}}\left( x \right)$ невозможно. Достаточно хорошее приближение позволяет получить КЭМ, основанный на последовательном уточнении параметров $v$ многомерной плотности распределения $q\left( {x,~v} \right),$ и процедура определения результирующей величины представляется в виде цепи Маркова с выбором параметров ВПР на каждом шаге. Для оценки параметров v данный метод использует результаты промежуточных статистических испытаний. Вектор параметров v определяется путем минимизации кросс-энтропии (KL-дивергенции) [28].

KL-дивергенция определяет меру близости двух произвольных ПР ${{p}_{1}}\left( x \right);$ ${{p}_{2}}\left( x \right){\text{:}}$

В данной работе в качестве ${{p}_{1}}\left( x \right);$ ${{p}_{2}}\left( x \right)$ принимаются соответственно оптимальная ПР $q{\text{*}}\left( x \right)~и~$ ее текущая оценка $q\left( {x;v} \right).$ KL-дивергенция ${{\mathbb{D}}_{{KL}}}\left( {q{\text{*}}(x),q\left( {x;v} \right)} \right)$ определяет близость между этими функциями. В предлагаемой работе функция $q\left( {x;v} \right)$ представляется плотностью многомерного нормального распределения с математическим ожиданием ${{\mu }_{q}}$ и матрицей корреляционных моментов ${{{{\Sigma }}}_{q}}.$ В такой постановке степень оптимальности, по существу, определяет вектор параметров $v = [{{\mu }_{q}};{{{{\Sigma }}}_{q}}]$). При этом процедура оптимизации может быть представлена в виде:

Подставляя для определения v* наилучшую ПР (5) в полученное выражение и рассматривая ${{\mathbb{P}}_{H}}$ как константу, получаем:

При использовании ВПР $q\left( x \right)$

где $W\left( x \right) = ~{{f\left( x \right)} \mathord{\left/ {\vphantom {{f\left( x \right)} {q\left( x \right)}}} \right. \kern-0em} {q\left( x \right)}}.$

КЭМ решает проблему оптимизации итеративно путем определения ряда промежуточных плотностей распределения $\{ q(x;{{v}_{k}}),k = 1,....,NT\} ,$ которые, как показано на рис. 1, постепенно приближаются к целевой плотности $q{\text{*}}\left( x \right),$ представляющей область существования редкого события.

На шаге $k$ оптимальная ПР $q{\text{*}}\left( x \right)$ может быть представлена оценкой $q\left( {x;v_{{k - 1}}^{*}} \right)$ с полученными на предыдущем шаге оптимальными параметрами $v_{{k - 1}}^{*}.$ При этом $~W\left( {x;v_{{k - 1}}^{*}} \right) = {{f(x)} \mathord{\left/ {\vphantom {{f(x)} {q\left( {x;v_{{k - 1}}^{*}} \right)}}} \right. \kern-0em} {q\left( {x;v_{{k - 1}}^{*}} \right)}}.$ Область ${{H}_{k}}$ промежуточного множества состояний системы определяется порогом ${{b}_{k}}{\text{:}}$

Порог ${{b}_{k}}$ вычисляется как θ-квантиль (например, дециль полученных в процессе статистических испытаний и отсортированных от наименьшего к наибольшему значениям пороговой функции ${{\psi }_{i}}\left( {{{x}_{i}}} \right).$ При этом моделирование выполняется согласно плотности распределения $q(x;{{v}_{{k - 1}}})$ с параметрами ${{v}_{{k - 1}}}.$

Полученный в результате решения задачи оптимизации на шаге $k - 1$ вектор $v{\text{*}}$ рассматривается как новое значение вектора параметров оптимальной ПР $q{\text{*}}\left( x \right),$ ${{v}_{k}} = v{\text{*}}.$ Начиная от исходного вектора параметров ${{v}_{0}},$ каждый последующий вектор ${{v}_{k}}$ определяется решением оптимизационной задачи (7), (9), приближаясь к оптимальной плотности распределения $q{\text{*}}\left( x \right) = \mathop {\lim }\limits_{k \to \infty } \left( {q(x;{{v}_{k}}){\text{\;}}} \right),$ которая, в свою очередь, является наилучшей оценкой оптимальной в кросс-энтропийном методе вспомогательной ПР.

При этом $W\left( {{{x}_{i}};{{v}_{{k - 1}}}} \right) = {{f(x)} \mathord{\left/ {\vphantom {{f(x)} {q\left( {x;{{v}_{{k - 1}}}} \right)}}} \right. \kern-0em} {q\left( {x;{{v}_{{k - 1}}}} \right)}}.$

Процедура повторяется до тех пор, пока ${{b}_{k}}$ не станет отрицательной, или, по крайней мере, $\theta \in [0.01,~0.1]$ испытаний не разместятся в искомой области редких событий [28]. Если число шагов до достижения критерия окончания итерационного процесса составляет m, то результирующая вероятность редкого события:

Как было упомянуто, вектор ${{v}_{0}}$ принимается равным исходным параметрам вероятностных распределений анализируемых случайных величин, т.е. при заданном $N$ и нулевом пороговом значении определяемая соотношением (8) область ${{H}_{0}},$ может и не содержать ни единого редкого события, что приводит к вырожденности оптимизационной процедуры (9). Именно поэтому область ${{H}_{k}}$ расширяется за счет ненулевого порогового значения ${{b}_{k}}.$ В результате в число событий-отказов войдут “почти отказы”. Отсюда отказы становятся менее редкими и идентифицируются в выборке.

Алгоритм $ECEM$

1. Ввод исходных данных: ${{v}_{0}},$ $\mathbb{C}{{\mathbb{V}}_{{target}}}.$ Начальные присвоения $k{\text{\;}} = {\text{\;}}0,$ ${{\delta }_{k}}$ = ∞.

2. k = k + 1. Генерация N реализаций состояний системы x на базе Гауссовского распределения с плотностью $q\left( {x,{{v}_{{k - 1}}}} \right).$

3. Формирование векторов: $J\left( {x,{{\delta }_{{k - 1}}}} \right);$ $q\left( {x,{{v}_{{k - 1}}}} \right);$ $\psi \left( x \right),$ а также

4. Определение ${{\delta }_{k}}$ при оптимизации: ${{\delta }_{k}} = \mathop {{\text{argmin}}}\limits_{{{{{\delta }}}_{k}} \in (0,{{{{\delta }}}_{{k - 1}}})} \left| {\mathbb{C}{{\mathbb{V}}_{{{{W}_{k}}}}} - \mathbb{C}{{\mathbb{V}}_{{target}}}} \right|.$

5. Определение ${{{v}}_{k}}$ при оптимизации: ${{{v}}_{k}} = \mathop {{\text{argmax}}}\limits_{{{{v}}_{k}}} \left[ {\frac{1}{N}\sum\nolimits_{i = 1}^N {{{W}_{k}}\left( {{{x}_{i}};{{{v}}_{{k - 1}}}} \right)\ln (q\left( {{{x}_{i}};{{{v}}_{k}}} \right))} } \right].$

6. Пп. 2–5 выполняются до тех пор, пока коэффициенты вариации ℂ𝕍 весов больше заданных $\mathbb{C}{{\mathbb{V}}_{{target}}}.$

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

ПРОВЕРОЧНЫЕ РАСЧЕТЫ

Для сравнения описанных процедур были выполнены расчеты показателей балансовой надежности пяти-узловой электрической схемы. Результаты расчетов представлены в табл. 1, где обозначены методы: ММК – классический метод Монте-Карло (при 10⁷ испытаниях); Моноцентр – метод подпространств с одной точкой формирования промежуточного подпространства; Полицентр – в качестве опорных точек формирования промежуточного подпространства ${{H}_{k}}$ принимается множество ${{H}_{{k - 1}}};$ SubSet – метод вложенных пространств с изменяющимися параметрами вероятностных распределений при формировании промежуточного подпространства; СЕ – классический кросс-энтропийный метод; ECE – модифицированный кросс-энтропийный метод.

Во втором столбце представлена вычисленная вероятность суммарного дефицита мощности ЭЭС. Следует отметить, что все представленные методы показывают приемлемую для практических расчетов точность редкого события – максимальное отклонение 3.3% (столбец Pr_Eps, %, метод СЕ) при точном значении вероятности Pr = = 2.23 × 10^–5 можно считать незначительным для малых вероятностей. Для сравнения можно упомянуть, что при таких вероятностях широко распространенная замена биномиального распределения распределением Пуассона имеет гораздо большую погрешность [29].

Последующие 3 столбца относятся к МО дефицита мощности. Здесь разброс (коэффициент вариации, столбец mD, Сv) результирующих величин значительно больше. Следует отметить, что относительная ошибка расчета МО дефицита мощности (коэффициент вариации, столбец mD, Сv)) превышает 25%. Увеличение точности возможно за счет увеличения объема промежуточной выборки. Однако это приводит к увеличению длительности расчетов (столбец t), которая в представленных расчетных процедурах зависит в основном от емких по времени преобразований вероятностных распределений (равномерное – Гауссовское – индивидуальное (маржинальное)). Если принять, что все случайные величины описываются одним и тем же нормальным распределением (с отличающимися параметрами), то это позволяет снизить длительность расчетов в несколько раз.

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

ЗАКЛЮЧЕНИЕ

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

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

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

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

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