Development and comparative analysis of mathematical models for the functioning of the regional power system of the Samara region

Full Text

Введение

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

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

1. Анализ известного математического описания функционирования региональной энергосистемы

В известных подходах к математическому моделированию функционирования энергетических систем в качестве исследовательского аппарата широкое распространение получили мультипликативно-степенные модели в виде производственных функций [9, 10]. В частности, в работах [1–3, 5, 7, 11], посвященных синтезу математических моделей региональной энергосистемы как многомерных производственных функций, используется нелинейная трехфакторная степенная зависимость вида
\[ \begin{equation}
Y_s ( t )=AK ^\alpha( t ) L ^\beta ( t )B ^\gamma ( t ),
\end{equation} \tag{1} \]
где $Y_s ( t )$ — выпуск электрической и тепловой энергии; $K ( t )$, $L ( t )$ и $B ( t )$ — капитальные, трудовые и топливные ресурсы соответственно.

Оценка параметров этой зависимости, прогноз и вычисление доверительных интервалов для результатов расчета в этих работах проводились на основе среднеквадратического оценивания коэффициентов линейной регрессионной модели
\[ \begin{equation}
\ln y_k =\ln A+\alpha \ln x_{1k} +\beta \ln x_{2k} +\gamma \ln x_{2k} +\eta_k,
\quad
\| \hat {\eta } \|^2=\sum\limits_{k=0}^{31} \hat {\eta }_k^2 \to \min,
\end{equation} \tag{2} \]
аппроксимирующей нелинейную зависимость (1), где $y_k$ — статистические данные по выпуску энергии; $x_{1k}$, $x_{2k} $ и $x_{3k} $ — статистические данные по динамике капитальных, трудовых и топливных ресурсов соответственно за период наблюдений с 1990 по 2021 годы в Самарской области.

Применяемые в этих работах методы математического моделирования имеют ряд существенных недостатков, негативно влияющих как на достоверность оценок основных показателей эффективности функционирования региональной энергосистемы, так и на точность прогноза, сделанного на основе построенной модели. Во-первых, недостаточно обоснован выбор вида и формы математической модели — линейной регрессионной трехфакторной модели (2), построенной на основе аппроксимации мультипликативно-степенной производственной функции (1). Для такой функциональной зависимости факторные эластичности
\[ \begin{equation*}
E_K =\cfrac{K ( t )}{Y_s ( t )} \cdot \cfrac{\partial Y_s }{\partial K}=\alpha,
\quad E_L =\frac{L ( t )}{Y_s ( t )} \cdot \cfrac{\partial Y_s }{\partial L}=\beta,
\quad
E_B =\cfrac{B (t )} {Y_s ( t )} \cdot \cfrac{\partial Y_s }{\partial B}=\gamma
\end{equation*} \]
являются параметрами модели, величинами априори постоянными, не изменяющимися на промежутке наблюдений, равном 31 году, что вызывает сомнение и требует подтверждения на основе статистических методов обработки результатов наблюдений [11].

Применение методов прикладного регрессионного анализа в задаче оценки параметров мультипликативно-степенной производственной функции (1) основано на минимизации среднеквадратичного отклонения результатов вычислений по модели $\hat {y}_k$ от результатов наблюдений $y_k$ в точках $t_k$ [11–14]:
\[ \begin{equation}
\sum\limits_{k=0}^{31} (y_k -\hat {y}_k )^2 = \| y-\hat {y} \|^2= \| e \|^2\to \min .
\end{equation} \tag{3} \]

Использование линейной аппроксимации на основе логарифмирования трехфакторной степенной зависимости (1) в регрессии (2) и критерий среднеквадратичного оценивания $\| \hat {\eta } \|^2=\sum _{k=0}^{31} \hat {\eta }_k^2 \to \min$ вносят смещение в оценки параметров $\alpha$, $\beta $ и $\gamma $, величина которого может оказаться недопустимо большой. В работах [1–8] эта проблема не исследуется.

Во-вторых, статистические данные, публикуемые в ежегодной отчетности региональных министерств и энергетических компаний, на основе которых строится математическая модель динамики выпуска продукции энергосистемы, представляют собой временной ряд наблюдений $y_k$, $k=0$, 1, 2, $\ldots$, 31, за период времени с 1990 по 2021 годы. Следовательно, для описания исследуемой динамики целесообразно выбирать модель из класса ковариационно-стационарных моделей временных рядов в форме стохастических разностных уравнений [16–20].

В-третьих, в работах [1–7] отсутствуют в необходимом объеме элементы статистического анализа результатов построения математической модели (1), а также выводы, сделанные на их основе. В частности, среднеквадратичное оценивание на основе минимизации суммы квадратов остатков $\| y-\hat{y} \|^2 =\| e \|^2\to \min$ для регрессионной модели дает наилучшие оценки только при выполнении условий теоремы Гаусса–Маркова [12, 14]. Если эти условия не выполняются, например, имеет место сильная корреляция между элементами вектора случайной помехи $\varepsilon_k$ и $\varepsilon_{k\pm p}$, $p\ne 0$, то применение метода наименьших квадратов (МНК) приводит к существенному занижению результатов расчета дисперсий и ковариаций среднеквадратичных оценок коэффициентов модели по сравнению с действительными значениями. Кроме того, оценки $\sigma_\varepsilon ^2 $ по остаточной сумме квадратов оказываются смещенными [12]. При этом создается ложное впечатление не только о точности оценок, но и о целесообразности проведения статистического анализа результатов оценивания, так как нет уверенности в корректности уровней значимости $t$- и $F$-критериев. Поэтому анализ остатков, выявление в них корреляции, например, на основе критерия Дарбина–Уотсона, является важнейшим элементом статистического анализа при построении математических моделей [13, 21]. Однако в работах [1–6] результаты таких исследований не приведены, что является серьезным недостатком.

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

В-четвертых, при использовании методов математического моделирования в опубликованных работах [1–6], посвященных системному анализу функционирования региональной энергосистемы, отсутствуют какие-либо исследования, связанные с проблемой мультиколлинеарности в регрессионном анализе [12, 14]. Игнорирование этой проблемы может привести к неустойчивости среднеквадратичных оценок, увеличению их дисперсий и, как следствие, к резкому снижению точности результатов расчета параметров модели и точности предсказания по модели.

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

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

2. Постановка задачи исследования и методы ее решения

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

При описании динамики выпуска продукции региональной энергосистемы может быть использована ковариационно-стационарная модель временного ряда в форме стохастического разностного уравнения в виде [20]
\[ \begin{equation}
y_k =\sum\limits_{j=1}^p \lambda _j y_{k-j} +
\sum\limits_{j=p+1}^{m+p+1} \lambda _j t_k^{j-( p+1 ) } + u (x_{1k} ,x_{2k} ,x_{3k} ,\bar {a} ) + \varepsilon _k,
\end{equation} \tag{4} \]
где $y_k =\sum _{j=1}^p \lambda _j y_{k-j} + \varepsilon _k$ — стохастическое разностное уравнение в форме модели авторегрессии порядка $p$; $P_m (t_k ) = \sum_{j=p+1}^{m+p+1} \lambda_j t_k^{j-(p+1)} $ — детерминированный тренд в форме многочлена степени $m$ относительно временной переменной $t_k$; $u (x_{1k},x_{2k},x_{3k},\bar {a} )$ — детерминированная трехфакторная функция с вектором параметров $\bar {a} = ( \lambda _{m+p+2} ,\lambda _{m+p+3} ,\ldots, \lambda _n )^\top$, $\varepsilon _k $ — последовательность независимых случайных величин с нулевым математическим ожиданием и одинаковыми дисперсиями $\sigma_\varepsilon^2$:
\[ \begin{equation*}
M [\varepsilon _k ]=0,
\quad
\operatorname{cov} [\varepsilon_k,\varepsilon _{k\pm p} ]=
\begin{cases}
0, & p \ne 0; \\
\sigma_\varepsilon^2, & p=0.
\end{cases}
\end{equation*} \]

При нелинейной зависимости $u ( x_{1k} ,x_{2k} ,x_{3k} ,\bar {a} )$ на основе ее линеаризации в окрестности точки $\bar {a}^{( i )}=
(\hat {\lambda }_{m+p+2}^{( i )},\hat {\lambda }_{m+p+3}^{( i )},\ldots,\hat {\lambda }_n^{( i )} )^\top$:
\[ \begin{equation*}
u (x_{1k} ,x_{2k} ,x_{3k} ,\bar {a} ) \approx
u_k^{( i )} + \sum\limits_{j=m+p+2}^n \frac{\partial u_k^{( i )}
}{\partial \lambda _j } (\lambda _j -\hat {\lambda }_j^{( i )} ),
\end{equation*} \]
где $u_k^{( i )} = ( x_{1k} ,x_{2k} ,x_{3k} ,\hat {\lambda }_{m+p+2}^{( i )} ,\hat {\lambda }_{m+p+3}^{( i )} ,\ldots,\hat {\lambda }_n^{( i )} )^\top$, ковариационно-стационарную модель временного ряда (4) можно привести к линейному виду $z_k^{( i )} = \lambda^\top f_k^{( i )} +\varepsilon _k$ или
\[ \begin{equation}
z^{( i )}=F^{( i )} \lambda +\varepsilon,
\end{equation} \tag{5} \]
где
\[ \begin{equation*}
z_k^{( i )} =y_k -u_k^{(i)} +\sum _{j=m+p+2}^n \frac{\partial u_k^{(i)} }{\partial
\lambda _j } \hat {\lambda }_j^{(i)}
\end{equation*} \]
— $k$-тый элемент вектора $z^{(i)}$;
\[ \begin{equation*}
f_k^{(i)} = \Bigl( y_{k-1} ,y_{k-2} ,\ldots, y_{k-p}, 1, t_k, t_k^2,\ldots, t_k^m, \cfrac{\partial u_k^{(i)} }{\partial \lambda _{m+p+2}},\ldots, \cfrac{\partial u_k^{(i)}}{\partial \lambda _n} \Bigr)
\end{equation*} \]
— $k$-тая строка матрицы $F^{(i)}$;
\[ \begin{equation*}
\lambda = (\lambda _1,\lambda _2 ,\ldots, \lambda_p, \lambda _{p+1}, \ldots, \lambda _{p+m+1}, \lambda _{p+m+2} ,\ldots, \lambda _n )^\top
\end{equation*} \]
— вектор искомых коэффициентов линейной регрессионной модели (5) при $k=p$, $p+1$, $\ldots$, $N-1$, где $N$ — объем выборки результатов наблюдений.

С учетом критерия среднеквадратичного оценивания (3) новое уточненное значение вектора $\hat {\lambda }^{( i+1 )}$ оценок коэффициентов разностного уравнения (4) находится по формуле
\[ \begin{equation}
\hat {\lambda }^{(i+1)}= [F^{( i )\top} F^{(i)} ]^{-1}F^{(i)\top}z^{( i )},\quad i=0, 1, 2, \ldots .
\end{equation} \tag{6} \]

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

К первой группе относятся наиболее простые трехфакторные регрессионные модели, образованные на основе производственных функций: линейной и мультипликативно-степенной:

• трехфакторная линейная регрессионная модель, построенная на основе линейной производственной функции (модель 1):
  \[ \begin{equation}
  y_k =\lambda _1 x_{1k} +\lambda _2 x_{2k} +\lambda _3 x_{3k} +\varepsilon_k,
  \quad k=0, 1, 2, \ldots, N-1;
  \end{equation} \tag{7} \]
• трехфакторная нелинейная регрессионная модель, построенная на основе нелинейной мультипликативно-степенной производственной функции (модель 2):
  \[ \begin{equation}
  y_k =\lambda _1 x_{1k}^{\lambda _2 } x_{2k}^{\lambda _3 } x_{3k}^{\lambda _4 } +\varepsilon _k,
  \quad k=0, 1, 2, \ldots, N-1.
  \end{equation} \tag{8} \]

Вторую группу образуют ковариационно-стационарные модели (4) в форме стохастических разностных уравнений временных рядов, включающие детерминированный полиномиальный тренд $P_m (t_k )$ и, в общем случае, нелинейную функцию входных воздействий (факторов) $u (x_{1k} ,x_{2k} ,x_{3k} \ )$:

• ковариационно-стационарная стохастическая модель в форме разностного уравнения с полиномиальным трендом второго порядка (модель 3):
  \[ \begin{equation}
  y_k =\lambda _1 y_{k-1} +\lambda _2 y_{k-2} +\lambda _3 +\lambda _4 t_k+\lambda _5 t_k^2 +\varepsilon _k,
  \quad
  k= 2, 3, 4, \dots, N-1;
  \end{equation} \tag{9} \]
• ковариационно-стационарная стохастическая модель в форме разностного уравнения с полиномиальным трендом нулевого порядка и линейной функцией входных воздействий/факторов (модель 4):
  \[ \begin{equation}
  y_k =\lambda _1 y_{k-1} +\lambda _2 +\lambda _3 x_{1k} +\lambda _4 x_{2k}+\lambda _5 x_{3k} +\varepsilon_k,
  \quad k= 1, 2, 3, \ldots, N-1;
  \end{equation} \tag{10} \]
• ковариационно-стационарная стохастическая модель в форме разностного уравнения с нелинейным трендом в форме мультипликативно-степенной трехфакторной функции (модель 5):
  \[ \begin{equation}
  y_k =\lambda _1 y_{k-1} +\lambda _2 y_{k-2} +\lambda _3 x_{1k}^{\lambda _4 } x_{2k}^{\lambda _5 } x_{3k}^{\lambda _6 } +\varepsilon _k ,
  \quad
  k= 2, 3, 4, \dots, N-1.
  \end{equation} \tag{11} \]

Третью группу образуют ковариационно-стационарные стохастические модели временных рядов, построенные на основе регрессионных моделей вида $y_k =u (x_{1k} ,x_{2k} ,x_{3k} ) + \eta_k$, в которых $u (x_{1k},x_{2k},x_{3k} )$ — производственная функция (линейная или мультипликативно-степенная) [10], а случайное возмущение $\eta_k$ описывается авторегрессионной моделью временного ряда первого или второго порядков соответственно:
\[ \begin{equation*}
\eta_k=\lambda _1 \eta _{k-1} +\varepsilon _k , \quad \eta _k =\lambda _1
\eta _{k-1} +\lambda _2 \eta _{k-2} +\varepsilon _k.
\end{equation*} \]

3. Разработка ковариационно-стационарных моделей временных рядов на основе линейной и нелинейной регрессии

При описании случайного возмущения моделью авторегрессии первого порядка имеем систему уравнений
\[ \begin{equation}
\begin{cases}
y_k =u ( x_{1k} ,x_{2k} ,x_{3k} )+\eta_k, & k=1,2,3,\ldots,N-1; \\
\eta _k =\lambda _1 \eta _{k-1} +\varepsilon _k,\quad |\lambda _1 |<1, & k=1,2,3,\ldots,N-1.
\end{cases}
\end{equation} \tag{12} \]

Второе уравнение в (12) можно доопределить для $\eta _0$. Найдем дисперсию $\sigma _\eta^2$. С учетом стационарности процесса и очевидных условий $M [\eta _k ]=0$ и $M [\eta _{k-1} \varepsilon _k ]=0$ имеем
\[ \begin{multline*}
\sigma_\eta^2 =\sigma ^2 [\eta _k ]=M\bigl[ \bigl(\eta_k -M [\eta _k ] \bigr)^2 \bigr]=
M [ \eta _k^2 ]=
M \bigl[ ( \lambda _1 \eta _{k-1} +\varepsilon _k )^2 \bigr]= {}
\\
{} =\lambda _1^2 M [\eta _{k-1}^2 ] + 2\lambda _1 M [ \eta_{k-1} \varepsilon _k ] +
M [\varepsilon _k^2 ] = \lambda_1^2 \sigma _\eta^2 +\sigma _\varepsilon^2.
\end{multline*} \]
Отсюда $\sigma_\eta ^2 = {\sigma _\varepsilon ^2 }/({1-\lambda _1^2 })$. Величину $\eta _0 $ будем искать из равенства $\mu \eta _0 =\varepsilon _0$ с учетом формулы $\sigma^2 [\mu \eta _0 ]=\mu^2 \sigma^2 [\eta _0 ] = \sigma^2 [\varepsilon _0 ]$. Отсюда $\mu ^2= {\sigma _\varepsilon^2}/{\sigma_\eta^2}=1-\lambda_1^2$, $\sqrt {1-\lambda _1^2} \eta _0 =\varepsilon _0$.

Тогда из соотношений
\[ \begin{equation*}
\begin{cases}
\sqrt {1-\lambda _1^2} \eta _0 =\varepsilon _0; \\
y_0 -u (x_{10} ,x_{20} ,x_{30} ) = \eta_0, & k=0,
\end{cases}
\; \;
\begin{cases}
\eta_k - \lambda_1 \eta_{k-1} =\varepsilon _k ; \\
\eta _k =y_k -u (x_{1k} ,x_{2k} ,x_{3k} ); \\
\eta _{k-1} =y_{k-1} -u (x_{1,k-1} ,x_{2,k-1} ,x_{3,k-1} ), \\
\qquad\qquad \qquad k=1, 2, 3,\ldots,N-1,
\end{cases}
\end{equation*} \]
получаем систему уравнений, описывающую временной ряд наблюдений:
\[ \begin{equation}
\begin{cases}
y_0 = y_0 -\sqrt {1-\lambda _1^2 } [ y_0 -u (x_{10},x_{20},x_{30} ) ]+\varepsilon _0; &
\\
y_k =u (x_{1k} ,x_{2k} ,x_{3k} ) + {} & \\
\qquad {}+ \lambda _1 [y_{k-1} - u (x_{1,k-1} ,x_{2,k-1} ,x_{3,k-1} ) ] + \varepsilon_k, &
k=1, 2, 3, \ldots, N-1.
\end{cases}
\end{equation} \tag{13} \]

Аналогично можно сформировать систему уравнений при описании случайного возмущения моделью авторегрессии второго порядка:
\[ \begin{equation}
\begin{cases}
y_k =u (x_{1k} ,x_{2k} ,x_{3k} )+\eta _k, & k=0, 1, 2, \ldots, N-1; \\
\eta_k = \lambda_1 \eta_{k-1} + \lambda_2 \eta_{k-2} + \varepsilon_k,
\quad |\lambda _2 |<1, & k=2, 3, 4, \ldots, N-1.
\end{cases}
\end{equation} \tag{14} \]

Найдем дисперсию $\sigma _\eta ^2 $ и ковариацию $\operatorname{cov} [\eta _k , \eta _{k-1} ]$. С учетом стационарности процесса и очевидных условий $M [\eta _k ]=0$, $M [\eta _{k-1} \varepsilon _k ]=0$ и $M [ \eta _{k-2} \varepsilon _k ]=0$ имеем
\[ \begin{multline*}
\operatorname{cov} [\eta _{k-1} ,\eta _{k-2} ] =
M [ \eta _{k-1} \eta_{k-2} ] =
M [ (\lambda _1 \eta _{k-2} +\lambda _2 \eta _{k-3} +\varepsilon _{k-1} ) \eta _{k-2} ]= {}
\\
{} =\lambda _1 \sigma _\eta ^2 +\lambda _2 M [ \eta _{k-2} \eta _{k-3} ] +
M [ \eta _{k-2} \varepsilon _{k-1} ] =
\lambda _1 \sigma _\eta^2 + \lambda_2 \operatorname{ cov} [\eta _{k-2} ,\eta _{k-3} ].
\end{multline*} \]
Отсюда $\operatorname{cov} [\eta _{k-1} ,\eta _{k-2} ] = {\lambda_1} \sigma_\eta^2/ ( {1-\lambda _2})$. Тогда
\[ \begin{multline*}
\sigma_\eta^2 =M [\eta_k^2 ] = M [ (\lambda _1 \eta _{k-1} +\lambda _2 \eta _{k-2} +\varepsilon _k )^2 ]={}
\\
{}
=\lambda _1^2 \sigma _\eta ^2 +\lambda _2^2 \sigma _\eta ^2 +\sigma_\varepsilon ^2 +2\lambda _1 \lambda _2
M[\eta _{k-1} \eta _{k-2} ]+ {}
\\
{}
+ 2\lambda _1 M [\eta _{k-1} \varepsilon _k ] + 2\lambda _2
M [ \eta _{k-2} \varepsilon _k ] =
\lambda _1^2 \sigma _\eta ^2 + \lambda _2^2 \sigma _\eta ^2 +\sigma _\varepsilon ^2 +\cfrac{2\lambda _1^2
\lambda _2 }{1-\lambda _2 }\sigma _\eta ^2 .
\end{multline*} \]
Отсюда после простых алгебраических преобразований получаем
\[ \begin{equation*}
\sigma _\eta ^2 =\cfrac{ (1-\lambda _2) \sigma _\varepsilon^2 }
{(1+\lambda _2 ) [ (1-\lambda _2 )^2 - \lambda _1^2 ]}.
\end{equation*} \]

Доопределим систему уравнений (14) для $\eta _0 $ и $\eta _1 $, которые будем искать из соотношений
\[ \begin{equation}
\begin{cases}
\mu_{11} \eta_0 = \varepsilon _0; \\
\mu_{21} \eta_0 + \mu _{22} \eta _1 =\varepsilon _1.
\end{cases}
\end{equation} \tag{15} \]
Из первого уравнения системы (15) получаем $\sigma^2 [\mu _{11} \eta _0 ] = \mu_{11}^2 \sigma^2 [\eta _0 ] = \sigma^2 [\varepsilon _0 ]$. Отсюда $\mu _{11}^2 = {\sigma _\varepsilon ^2 }/{\sigma _\eta ^2 }=( 1+\lambda _2 ) [ ( 1-\lambda _2 )^2 - \lambda _1^2 ] /(1-\lambda _2)$, следовательно,
\[ \begin{equation}
\mu _{11} =\sqrt { ( 1+\lambda _2 ) [( 1-\lambda _2 )^2-\lambda _1^2 ]/(1-\lambda _2)}.
\end{equation} \tag{16} \]

Элементы $\mu _{21}$ и $\mu _{22}$ найдем из системы равенств
\[ \begin{equation*}
\begin{cases}
\sigma^2 [\mu _{21} \eta _0 +\mu _{22} \eta _1 ] = \sigma^2 [\varepsilon _1 ];
\\
\operatorname{cov} [ \varepsilon _0 , \varepsilon _1 ]=0.
\end{cases}
\end{equation*} \]
Из первого равенства получаем $\mu _{21}^2 \sigma _\eta ^2 +2\mu _{21} \mu_{22} \operatorname{ cov} [\eta _0 ,\eta _1 ] + \mu _{22}^2 \sigma _\eta ^2 =\sigma _\varepsilon^2$. Из второго — $\operatorname{cov} [ \varepsilon _0 ,\varepsilon _1 ] = \operatorname{cov} [\mu _{11} \eta_0, \mu _{21} \eta _0 +\mu _{22} \eta _1 ] = \mu _{11} \mu _{21} \sigma _\eta ^2 +\mu _{11} \mu _{22}{ \operatorname{cov} [ \eta _0 ,\eta _1 ] = 0}$. Отсюда имеем
\[ \begin{equation*}
\mu _{21} =-\cfrac{\operatorname{cov} [ \eta _0 ,\eta _1 ]}{\sigma _\eta ^2}\mu _{22} =
-\cfrac{\lambda _1 }{1-\lambda _2 }\mu _{22},
\end{equation*} \]
\[ \begin{equation*}
\cfrac{\lambda _1^2 }{ ( 1-\lambda _2 )^2} \mu _{22}^2 \sigma_\eta^2 -
\cfrac{2\lambda _1 }{1-\lambda _2 }\mu _{22}^2
\cfrac{\lambda _1 }{1-\lambda _2 }\sigma _\eta ^2 +
\mu _{22}^2 \sigma _\eta ^2 =\sigma_\varepsilon^2,
\end{equation*} \]
\[ \begin{equation*}
\Bigl[1-\cfrac{\lambda _1^2 }{ (1-\lambda _2 )^2} \Bigr] \mu_{22}^2 =
\cfrac{\sigma _\varepsilon ^2}{\sigma _\eta ^2} =
\cfrac{ (1+\lambda _2 ) [ (1-\lambda _2 )^2 -\lambda_1^2 ]}{1-\lambda _2},
\end{equation*} \]
\[ \begin{equation}
\mu _{22} =\sqrt {1-\lambda _2^2} ,
\quad
\mu _{21} =-\cfrac{\lambda _1 }{1-\lambda _2 }\mu _{22} =-\lambda _1
\cfrac{\sqrt {1+\lambda _2 } }{\sqrt {1-\lambda _2 }}.
\end{equation} \tag{17} \]
Тогда из соотношений
\[ \begin{equation*}
\begin{cases}
\mu_{11} \eta _0 =\varepsilon_0; \quad\qquad \quad
y_0 - u (x_{10} ,x_{20} ,x_{30} )=\eta_0, & k=0; \\
\mu _{21} \eta _0 +\mu _{22} \eta _1 =\varepsilon _1; \quad
y_1 - u (x_{11} ,x_{21} ,x_{31} )=\eta_1, & k=1; \\
\eta_k -\lambda _1 \eta _{k-1} -\lambda _2 \eta _{k-2} =\varepsilon_k; \;\;
\eta _k =y_k -u (x_{1k} ,x_{2k} ,x_{3k} ); & \\
\eta _{k-1} =y_{k-1} -u (x_{1,k-1} ,x_{2,k-1} ,x_{3,k-1} ); & \\
\eta _{k-2} =y_{k-2} -u (x_{1,k-2} ,x_{2,k-2} ,x_{3,k-2} ), & k=2, 3, \ldots,N-1,
\end{cases}
\end{equation*} \]
получаем систему уравнений, описывающую временной ряд наблюдений:
\[ \begin{equation}
\begin{cases}
y_0 =y_0 -\mu _{11} [ y_0 -u (x_{10} ,x_{20} ,x_{30} ) ]+\varepsilon _0;
\\
y_1 =y_1 -\mu _{21} [ y_0 -u (x_{10} ,x_{20} ,x_{30} ) ]-
\mu _{22} [ y_1 -u (x_{11} ,x_{21} ,x_{31} ) ]+\varepsilon _1; \\
y_k =u ( x_{1k} ,x_{2k} ,x_{3k} ) +
\lambda _1 [ y_{k-1} -u ( x_{1,k-1} ,x_{2,k-1} ,x_{3,k-1} ) ]+ {}\\
\qquad {}
+\lambda _2 [ y_{k-2} -u ( x_{1,k-2} ,x_{2,k-2} ,x_{3,k-2} ) ]+\varepsilon _k, \quad k=2, 3, \ldots, N-1,
\end{cases}
\end{equation} \tag{18} \]
в которой зависимости $\mu _{11} (\lambda _1 ,\lambda _2 )$, $\mu _{21} (\lambda _1 ,\lambda _2 )$ и $\mu _{22} (\lambda _2 )$ определяются формулами (16) и (17).

При линейной $u (x_1 ,x_2 ,x_3 )=\alpha x_1 +\beta x_2 +\gamma x_3$ и мультипликативно-степенной $u (x_1 ,x_2 ,x_3 ) = Ax_1^\alpha x_2^\beta x_3^\gamma $ производственных функциях системы уравнений (13) и (18) принимают следующий вид (модели 7–9):
\[ \begin{equation}
\begin{cases}
y_0 =y_0 -\sqrt {1-\lambda _1^2 } (y_0 -\lambda _2 x_{10} -\lambda_3 x_{20} -\lambda _4 x_{30} )+\varepsilon _0; \\
y_k =\lambda _2 x_{1k} +\lambda _3 x_{2k} +\lambda _4 x_{3k} + {}\\
{}
\qquad {} +\lambda _1 (y_{k-1} -\lambda _2 x_{1,k-1} -\lambda _3 x_{2,k-1} -\lambda _4 x_{3,k-1} )+\varepsilon _k,
\\
\hspace{5.6cm}
k=1, 2, 3, \ldots, N-1,
\end{cases}
\end{equation} \tag{19} \]
\[ \begin{equation}
\begin{cases}
y_0 =y_0 -\sqrt {1-\lambda _1^2 } (y_0 -\lambda _2 x_{10}^{\lambda_3 } x_{20}^{\lambda _4 } x_{30}^{\lambda _5 } )+\varepsilon _0;
\\
y_k =\lambda _2 x_{1k}^{\lambda _3 } x_{2k}^{\lambda _4 } x_{3k}^{\lambda_5 } +
\lambda _1 ( y_{k-1} -\lambda _2 x_{1,k-1}^{\lambda _3 } x_{2,k-1}^{\lambda _4 } x_{3,k-1}^{\lambda _5 } ) + \varepsilon _k, \\
\hspace{6cm}
k=1, 2, 3, \ldots, N-1,
\end{cases}
\end{equation} \tag{20} \]
\[ \begin{equation}
\begin{cases}
y_0 = y_0 -\mu _{11} (y_0 -\lambda _3 x_{10} -\lambda _4 x_{20} -\lambda _5 x_{30} )+\varepsilon _0; \\
y_1 = y_1 -\mu _{21} (y_0 -\lambda _3 x_{10} -\lambda _4 x_{20} -\lambda _5 x_{30} )- {}\\
\qquad {}-\mu _{22} (y_1 -\lambda _3 x_{11} -\lambda _4 x_{21} -\lambda _5 x_{31} ) + \varepsilon_1; \\
y_k =\lambda _3 x_{1k} +\lambda _4 x_{2k} +\lambda _5 x_{3k} + {}\\
\qquad {} +\lambda _1 ( y_{k-1} -\lambda _3 x_{1,k-1} -\lambda _4 x_{2,k-1} -\lambda _5 x_{3,k-1} )+ {}
\\
\qquad {} +\lambda _2 ( y_{k-2} -\lambda _3 x_{1,k-2} -\lambda _4 x_{2,k-2} - \lambda _5 x_{3,k-2} )+\varepsilon _k,
\\
\hspace{5.7cm}
k=2, 3, 4, \ldots, N-1,
\end{cases} \!
\end{equation} \tag{21} \]
\[ \begin{equation}
\begin{cases}
y_0 =y_0 -\mu _{11} (y_0 -\lambda _3 x_{10}^{\lambda _4 } x_{20}^{\lambda _5 } x_{30}^{\lambda _6 } )+\varepsilon_0;
\\
y_1 =y_1 -\mu _{21} ( y_0 -\lambda _3 x_{10}^{\lambda _4 } x_{20}^{\lambda _5 } x_{30}^{\lambda _6 } )- {}
\\
\qquad
{}-\mu _{22} (y_1-\lambda _3 x_{11}^{\lambda _4 } x_{21}^{\lambda _5 } x_{31}^{\lambda _6 } )+\varepsilon _1 ; \\
y_k =\lambda _3 x_{1k}^{\lambda _4 } x_{2k}^{\lambda _5 } x_{3k}^{\lambda_6} + \lambda _1
( y_{k-1} -\lambda _3 x_{1,k-1}^{\lambda _4 } x_{2,k-1}^{\lambda _5 } x_{3,k-1}^{\lambda _6 } )+ {}
\\
\qquad {}
+\lambda _2 (y_{k-2} -\lambda _3 x_{1,k-2}^{\lambda _4 } x_{2,k-2}^{\lambda _5 } x_{3,k-2}^{\lambda _6 } )+\varepsilon _k, \\
\hspace{5.7cm}
k=2, 3, 4, \ldots, N-1,
\end{cases}
\end{equation} \tag{22} \]
где зависимости $\mu _{11} (\lambda _1 ,\lambda _2 )$, $\mu _{21} (\lambda _1 ,\lambda _2 )$ и $\mu _{22} (\lambda_2 )$ также определяются формулами (16) и (17).

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

Модели (7), (9) и (10) — линейные, а модели (8), (11), (19)–(22) — нелинейные относительно параметров, подлежащих оценке на основе результатов наблюдений в форме статистических данных, публикуемых в ежегодной отчетности региональных министерств и энергетических компаний. Отметим также, что во всех представленных выше моделях (7)–(11), (19)–(22) предполагается, что случайные возмущения $\varepsilon _k$, $k=0, 1, 2, \ldots, N-1$, в результатах наблюдений имеют нулевое математическое ожидание и не коррелируют друг с другом:
\[ \begin{equation*}
M [\varepsilon _k ]=0,
\quad
M [\varepsilon _k \varepsilon _{k\pm s} ]=
\begin{cases}
\sigma_\varepsilon^2, & s=0; \\
0, & s\ne 0.
\end{cases}
\end{equation*} \]
Поэтому в качестве критерия среднеквадратичного оценивания параметров моделей используется минимизация суммы квадратов остатков (3). Проверка обоснованности выбора такого критерия строится на основе анализа остатков, выявления в них корреляции и статистики Дарбина–Уотсона [13, 21].

Применение методов регрессионного анализа позволяет находить, а для нелинейных моделей уточнять оценки коэффициентов математической модели по формуле (6). При этом элементы вектора $z$ и элементы $f_{k,j}$ вектора-строки $f_k$ матрицы $F$ описываются следующими соотношениями.

1. Для трехфакторной линейной регрессионной модели (7), построенной на основе линейной производственной функции:
   \[ \begin{equation*}
   \begin{cases}
   z_k =y_{k-1}; \\
   f_k = ( x_{1,k-1} ,x_{2,k-1} ,x_{3,k-1} ), & k=1, 2, 3, \ldots, N.
   \end{cases}
   \end{equation*} \]
2. Для трехфакторной линейной регрессионной модели (8), построенной на основе мультипликативно-степенной производственной функции:
   \[ \begin{equation*}
   \begin{cases}
   u_k^{( i )} =\lambda _1^{( i )} x_{1k}^{\lambda_2^{( i )} } x_{2k}^{\lambda _3^{( i )} }
   x_{3k}^{\lambda _4^{( i )} }, \qquad k=0, 1, 2, \ldots, N-1;
   \\
   z_k^{(i)} =y_{k-1} +u_{k-1}^{(i)} (\lambda_2^{(i)} \ln x_{1,k-1} +\lambda _3^{(i)} \ln
   x_{2,k-1} +\lambda _4^{(i)} \ln x_{3,k-1} );
   \\
   f_k^{(i)} =u_{k-1}^{(i)} ( ({\lambda _1^{(i)} })^{-1}, \ln x_{1,k-1} , \ln x_{2,k-1} , \ln
   x_{3,k-1} ), \quad
   k=1, 2, 3, \ldots, N.
   \end{cases}
   \end{equation*} \]
3. Для линейной ковариационно-стационарной модели (9) в форме разностного уравнения с полиномиальным трендом второго порядка:
   \[ \begin{equation*}
   \begin{cases}
   z_k =y_{k+1} ; \\
   f_k = (y_k ,y_{k-1} ,1,t_{k+1} ,t_{k+1}^2 ), &
   k=1, 2, 3, \ldots, N-2.
   \end{cases}
   \end{equation*} \]
4. Для ковариационно-стационарной модели (10) в форме разностного уравнения с полиномиальным трендом нулевого порядка и линейной функцией входных воздействий (факторов):
   \[ \begin{equation*}
   \begin{cases}
   z_k =y_k ; \\
   f_k = (y_{k-1} ,1,x_{1k} ,x_{2k} ,x_{3k} ), &
   k=1, 2, 3, \ldots, N-1.
   \end{cases}
   \end{equation*} \]
5. Для ковариационно-стационарной модели (11) в форме разностного уравнения с нелинейным трендом в форме мультипликативно-степенной трехфакторной функции:
   \[ \begin{equation*}
   \begin{cases}
   u_k^{(i)} = \lambda _3^{(i)} x_{1k}^{\lambda _4^{(i)} } x_{2k}^{\lambda _5^{(i)} }
   x_{3k}^{\lambda _6^{(i)} }, \qquad
   k = 2, 3, 4, \ldots, N-1;
   \\
   z_k^{(i)} = y_{k+1} +u_{k+1}^{(i)} (\lambda_4^{(i)} \ln x_{1,k+1} + \lambda _5^{(i)} \ln
   x_{2,k+1} + \lambda _6^{(i)} \ln x_{3,k+1} );
   \\
   f_k^{(i)} = ( y_k , y_{k-1} , u_{k+1}^{( i )} /\lambda _3^{(i)} , u_{k+1}^{(i)} \ln
   x_{1,k+1} ,u_{k+1}^{(i)} \ln x_{2,k+1} ,u_{k+1}^{( i )} \ln x_{3,k+1} ),
   \\
   \hspace{8.1cm}
   k = 1, 2, 3, \ldots, N-2.
   \end{cases}
   \end{equation*} \]
6. Для ковариационно-стационарной модели (19), построенной на основе линейной регрессионной модели со случайным возмущением, описываемым авторегрессией первого порядка:
   \[ \begin{equation*}
   \begin{cases}
   u_k^{(i)} =\lambda _2^{(i)} x_{1k} +\lambda_3^{(i)} x_{2k} +\lambda _4^{(i)} x_{3k}
   , \qquad k=0, 1, 2, \ldots, N-2; \\
   z_k^{(i)} = y_k -\lambda _1^{(i)} u_{k-1}^{(i)}; \\
   f_k^{(i)} = ( y_{k-1} -u_{k-1}^{(i)}, x_{1k} -\lambda _1^{(i)} x_{1,k-1} ,x_{2k} -\lambda_1^{(i)} x_{2,k-1} ,x_{3k} -\lambda _1^{(i)}x_{3,k-1} ),
   \\
   \hspace{8.1cm}
   k = 1, 2, 3, \ldots, N-1.
   \end{cases}
   \end{equation*} \]
7. Для ковариационно-стационарной модели (20), построенной на основе нелинейной регрессионной модели со случайным возмущением, описываемым авторегрессией первого порядка:
   \[ \begin{equation*}
   \begin{cases}
   u_k^{(i)} =\lambda _2^{(i)} x_{1k}^{\lambda _3^{(i)} } x_{2k}^{\lambda _4^{(i)} } x_{3k}^{\lambda _5^{(i)} }, \qquad k=0, 1, 2, \ldots, N-1;
   \\
   z_k^{(i)} =y_k +u_k^{(i)} (\lambda_3^{(i)} \ln x_{1k} +\lambda _4^{(i)} \ln x_{2k} +\lambda_5^{(i)} \ln x_{3k} )- {}
   \\
   \qquad {}
   -\lambda _1^{(i)} u_{k-1}^{(i)} (1+\lambda _3^{(i)} \ln x_{1,k-1} +\lambda _4^{( i )} \ln x_{2,k-1} +\lambda _5^{(i)} \ln x_{3,k-1} );
   \\
   f_k^{(i)}= \bigl( y_{k-1} -u_{k-1}^{(i)}, ( u_k^{( i )} -\lambda _1^{(i)} u_{k-1}^{(i)} )/ {\lambda _2^{(i)}}, \\
   \qquad\qquad
   u_k^{(i)} \ln x_{1k} -\lambda _1^{( i )} u_{k-1}^{(i)} \ln x_{1,k-1},
   u_k^{(i)} \ln x_{2k} - \lambda _1^{( i )} u_{k-1}^{(i)} \ln x_{2,k-1}, \\
   \qquad\qquad \qquad\qquad \qquad\qquad 
   u_k^{( i )} \ln x_{3k} - \lambda _1^{(i)} u_{k-1}^{( i )} \ln x_{3,k-1} \bigr),
   \\
   \hspace{8cm}
   k = 1, 2, 3, \ldots, N-1.
   \end{cases}
   \end{equation*} \]
8. Для ковариационно-стационарной модели (21), построенной на основе линейной регрессионной модели со случайным возмущением, описываемым авторегрессией второго порядка:
   \[ \begin{equation*}
   \begin{cases}
   u_k^{(i)} =\lambda _3^{(i)} x_{1k} +\lambda _4^{(i)} x_{2k} +\lambda _5^{(i)} x_{3k} , \qquad k=0, 1, 2,\ldots, N-2;
   \\
   z_k^{(i)} =y_{k+1} -\lambda _1^{(i)} u_k^{(i)} - \lambda _2^{(i)} u_{k-1}^{(i)};
   \\
   f_k^{(i)} = \bigl(y_k -u_k^{(i)}, y_{k-1} -u_{k-1}^{(i)}, x_{1,k+1} -\lambda_1^{(i)} x_{1,k} -\lambda _2^{(i)} x_{1,k-1} ,
   \\
   \qquad\qquad
   x_{2,k+1} -\lambda _1^{(i)} x_{2,k} - \lambda _2^{(i)} x_{2,k-1},
   x_{3,k+1} -\lambda_1^{(i)} x_{3,k} - \lambda _2^{(i)} x_{3,k-1}
   \bigr),
   \\
   \hspace{8cm}
   k = 1, 2, 3, \ldots, N-2.
   \end{cases}
   \end{equation*} \]
9. Для ковариационно-стационарной модели (22), построенной на основе нелинейной регрессионной модели со случайным возмущением, описываемым авторегрессией второго порядка:
   \[ \begin{equation*}
   \begin{cases}
   u_k^{(i)} =\lambda _3^{(i)} x_{1k}^{\lambda _4^{(i)} } x_{2k}^{\lambda _5^{(i)} } x_{3k}^{\lambda _6^{(i)}}, \qquad k=0, 1, 2, \ldots, N-1;
   \\
   z_k^{(i)} =y_{k+1} +u_{k+1}^{(i)} ( \lambda_4^{(i)} \ln x_{1,k+1} +\lambda _5^{(i)} \ln
   x_{2,k+1} +\lambda _6^{(i)} \ln x_{3,k+1} ) - {}
   \\
   \qquad {}
   -\lambda _1^{(i)} u_k^{(i)} ( 1+\lambda_4^{(i)} \ln x_{1,k} +\lambda _5^{(i)} \ln x_{2,k}
   +\lambda _6^{(i)} \ln x_{3,k} )- {}
   \\
   \qquad {}
   -\lambda _2^{(i)} u_{k-1}^{(i)} ( 1+\lambda _4^{(i)} \ln x_{1,k-1} +\lambda _5^{( i )} \ln x_{2,k-1} +\lambda _6^{(i)} \ln x_{3,k-1} );
   \\
   f_k^{(i)} =\bigl( y_k -u_k^{(i)}, y_{k-1} -u_{k-1}^{(i)}, ( u_{k+1}^{(i)} -\lambda_1^{(i)} u_k^{(i)} -\lambda _2^{(i)} u_{k-1}^{(i)} )/ \lambda_3^{(i)} ,
   \\
   \qquad\qquad
   u_{k+1}^{(i)} \ln x_{1,k+1} -\lambda _1^{(i)} u_k^{(i)} \ln x_{1,k} -\lambda _2^{(i)}
   u_{k-1}^{(i)} \ln x_{1,k-1},
   \\
   \qquad\qquad
   u_{k+1}^{(i)} \ln x_{2,k+1} -\lambda _1^{(i)} u_k^{(i)} \ln x_{2,k} -\lambda _2^{(i)} u_{k-1}^{(i)} \ln x_{2,k-1},
   \\
   \qquad\qquad \qquad\qquad
   u_{k+1}^{(i)} \ln x_{3,k+1} -\lambda _1^{( i )} u_k^{(i)} \ln x_{3,k} -\lambda _2^{(i)}
   u_{k-1}^{(i)} \ln x_{3,k-1} \bigr),
   \\
   \hspace{8cm}
   k = 1, 2, 3, \ldots, N-2.
   \end{cases}
   \end{equation*} \]

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

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

Таблица 1. Результаты расчета параметров математических моделей
[Results of the calculation of mathematical model parameters]

Model		Model coefficients						$Q_{res}$	$s$, %	cond	DW
nos.	Eq.	$\lambda_1$	$\lambda_2$	$\lambda_3$	$\lambda_4$	$\lambda_5$	$\lambda _6$
1	(7)	$-0.038$	$-0.048$	$1.103$				0.022	4.2	214	0.8
2	(8)	1.008	$-0.167$	$-0.058$	1.035			0.017	3.7	71	1.0
3	(9)	0.921	$-0.199$	108.2	$-0.1064$	$2.6\cdot 10^{-5}$		0.032	5.5	$1.5\cdot 10^7$	2.2
4	(10)	0.276	0.070	$-0.047$	$-0.053$	$-0.750$		0.011	3.2	4535	1.6
5	(11)	0.314	$-0.038$	0.730	$-0.158$	$-0.099$	1.071	0.011	3.3	2469	1.6
6	(19)	0.785	$-0.009$	0.009	0.953			0.013	3.4	170	1.9
7	(20)	0.684	0.965	$-0.128$	0.016	0.945		0.012	3.3	77	1.9
8	(21)	0.885	$-0.123$	0.0099	0.017	0.894		0.013	3.5	769	2.1
9	(22)	0.818	$-0.108$	0.900	$-0.063$	0.037	0.853	0.012	3.4	203	2.1

В последних четырех столбцах табл. 1 приведены следующие значения: $Q_{res} = \| y-\hat {y} \|^2$ — сумма квадратов отклонений результатов расчета $\hat{y}_k$ по модели от результатов наблюдений; $s= ( \| y-\hat {y} \| / \| y \| ) \cdot 100\,\%$ — оценка относительного отклонения модели от эксперимента; $\rm cond$ — число обусловленности матрицы нормальной системы уравнений $F^{(i)\top}F^{( i )}$, характеризующее устойчивость вычисления среднеквадратичных оценок параметров модели; $\mathrm{DW}=2(1-r)$ — статистика Дарбина–Уотсона, на основе которой можно сделать вывод о наличии автокорреляции первого порядка в случайном возмущении $\varepsilon_k$, где $r$ — выборочный коэффициент корреляции.

Статистический анализ полученных результатов, приведенных в табл. 1, позволяет сделать следующие выводы.

Во-первых, построенные модели 4–9 (см. формулы (10), (11), (19)–(22)) с достаточно высокой степенью адекватности аппроксимируют результаты наблюдений. Причем графики функций, описывающих динамику выпуска продукции энергопроизводств, построенные на основе моделей 4 и 5 (сплошная линия на рис. 1) и моделей 6–9 (штриховая линия на рис. 1), практически совпадают. Точки на рис. 1 отображают статистические данные выпуска продукции энергопроизводств в относительных единицах.

Рис. 1. Графики функций, описывающих динамику выпуска продукции энергопроизводств: точки — статистические данные выпуска продукции; сплошная линия — данные на основе моделей 4 и 5 (формулы (10) и (11)); штриховая линия — данные на основе моделей 6–9 (формулы (19)–(22))
[Figure. 1. Plots of functions describing the dynamics of production output in energy industries: points — statistical data on production output; solid line — data based on models 4 and 5 (Eqs. (10) и (11)); dashed line — data based on models 6–9 (Eqs. (19)–(22))]

Во-вторых, с учетом плохой обусловленности матрицы нормальной системы уравнений $F^{( i )\top}F^{( i )}$ при вычислении среднеквадратичных оценок параметров математической модели (9) в форме разностного уравнения с полиномиальным трендом второго порядка можно сделать вывод о наличии мультиколлинеарности, отрицательные последствия которой проявляются в виде неустойчивости как оценок, так и самой процедуры оценивания, а также в резком увеличении дисперсии оценок параметров и уменьшении точности предсказания по модели [12, 14]. Вследствие этого данная модель не может быть рекомендована к практическому применению при описании динамики выпуска продукции региональной энергосистемы и достоверного прогноза на ее основе.

В-третьих, из анализа остатков $e =y_k -\hat {y}_k$, $k=0,1,2,\ldots,31$, и оценки DW следует, что практически во всех моделях временного ряда наблюдений (из первой и второй групп по описанной выше систематизации) имеет место автокорреляция между элементами случайной величины $\varepsilon _k $. В связи с этим среднеквадратичное оценивание параметров модели на основе классического метода наименьших квадратов (МНК): $\| e \|^2= \| y-\hat {y} \|^2\to \min$ приводит не только к неэффективности этих оценок, но и к существенному занижению оценок дисперсий результатов расчета, в силу чего оценки дисперсии $\sigma _\varepsilon ^2$ оказываются смещенными и создается неверное впечатление о точности прогноза по модели [12, 14].

В-четвертых, статистический анализ построенных моделей (19)–(22) из третьей группы показал нецелесообразность использования при описании временного ряда случайных возмущений в результатах наблюдений модели авторегрессии второго порядка $\eta _k =\lambda _1 \eta _{k-1} +\lambda _2 \eta _{k-2} +\varepsilon _k$, а вполне достаточно ограничиться авторегрессией первого порядка $\eta _k =\lambda _1 \eta _{k-1} + \varepsilon _k$. То есть отдается предпочтение моделям (19) и (20).

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

5. Сравнительный анализ разработанных математических моделей на основе оценки погрешности прогноза

Сравнительный анализ исследуемых моделей был выполнен на основе данных, собранных за 20 лет ($N = 20$). Прогнозы были рассчитаны на сроки $p = 1$, 2, 3, 4 и 5 лет, начиная с 2010 года и до $2022-p$ года. Усредненные результаты погрешностей прогнозирования представлены в табл. 2. Анализ ошибок прогнозов при использовании моделей 1–9 позволяет сделать выбор в пользу модели 7, описываемой формулой (20).

Таблица 2. Погрешности прогноза выпуска продукции энергопроизводств, сделанного на основе исследуемых математических моделей (в процентах)
[Forecast errors in the production output of energy industries based on the studied mathematical models (in percent)]

Number of years for which the Model numbers forecast has been made	Model numbers
	1	2	3	4	5	6	7	8	9
1	7.6	4.9	6.8	4.4	4.5	4.9	3.7	5.8	5.4
2	9.1	7.3	9.8	7.2	7.0	7.4	5.5	7.9	7.9
3	10.3	9.6	11.6	10.1	9.1	9.5	7.2	10.6	9.8
4	9.9	11.0	15.5	12.7	11.4	11.1	9.2	11.3	11.0
5	9.0	12.1	15.2	15.1	13.4	12.1	11.2	12.1	12.5

Результаты сравнения точности прогноза на основе известной [1–3] мультипликативно-степенной модели (модель 2, формула (8)) и модели временного ряда в форме разностного уравнения (модель 7, формула (20)) представлены на рис. 2. Видно, что погрешность прогноза $\delta$, сделанного на основе модели 7, существенно меньше (в среднем на 27%) погрешности прогноза, сделанного на основе модели 2. Также видно, что из-за некорректности применения классического метода наименьших квадратов при оценке параметров модели 2 вследствие смещения оценки дисперсии случайного возмущения в результатах наблюдения величина оценки относительной предельной погрешности прогноза $\Delta$ (доверительный интервал) оказывается завышенной.

Рис. 2. Относительная ошибка прогноза $\delta$ (%) и относительная предельная ошибка прогноза $\Delta$ (%) объемов производства в отраслях энергетики, рассчитанных по моделям 2 и 7
[Figure. 2. Relative forecast error $\delta$ (%) and relative limit forecast error $\Delta$ (%) of production output in energy industries calculated by models 2 and 7]

При выборе формы трехфакторной функциональной зависимости (1), описывающей динамику выпуска продукции региональной энергосистемы, предполагалось, что факторные эластичности для наблюдаемого процесса есть величины постоянные, равные параметрам модели $\alpha$, $\beta$ и $\gamma$. Однако это предположение ни в одной из работ [1–8] не было обосновано.

Для устранения этого недостатка и обоснования независимости от времени факторных эластичностей построенной модели
\[ \begin{equation*}
y_k =\lambda _2 x_{1k}^{\lambda _3 } x_{2k}^{\lambda _4 } x_{3k}^{\lambda _5} +\lambda _1 ( y_{k-1} -\lambda _2 x_{1,k-1}^{\lambda _3 }
x_{2,k-1}^{\lambda _4 } x_{3,k-1}^{\lambda _5 } ) + \varepsilon _k
\end{equation*} \]
рассмотрим аппроксимацию зависимостей $\alpha ( t )$, $\beta ( t )$ и $\gamma ( t )$ многочленами нулевой степени $\alpha_0$, $\beta _0$ и $\gamma _0$ и проверим гипотезу об адекватности таких моделей результатам расчетов, представленных во втором, четвертом и шестом столбцах табл. 3 соответственно.

Таблица 3. Результаты расчетов коэффициентов эластичности и оценки точности этих расчетов
[Results of the calculations of elasticity coefficients and evaluation of the accuracy of these calculations]

Observation period	$\alpha$	$s [ \alpha ]$	$\beta$	$s [\beta]$	$\gamma$	$s [\gamma]$
1991–2011	0.093	0.198	0.061	0.110	0.881	0.258
1992–2012	0.159	0.181	0.105	0.110	0.658	0.298
1993–2013	0.155	0.181	0.102	0.110	0.681	0.258
1994–2014	0.210	0.176	0.142	0.113	0.562	0.238
1995–2015	0.188	0.136	0.118	0.095	0.548	0.225
1996–2016	0.226	0.109	0.106	0.073	0.457	0.207
1997–2017	0.103	0.136	0.107	0.075	0.499	0.189
1998–2018	$-0.072$	0.085	0.110	0.085	0.650	0.166
1999–2019	$-0.116$	0.059	0.122	0.065	0.700	0.143
2000–2020	$-0.109$	0.066	0.121	0.060	0.688	0.124
2001–2021	$-0.167$	0.064	0.195	0.077	0.673	0.120

Обозначим в общем виде факторную эластичность через $p\in \{\alpha,\beta,\gamma\}$ и рассмотрим регрессионную модель в виде многочлена нулевой степени
\[ \begin{equation}
p_k =c_0 +\eta _k, \quad k=\overline{1,n_p},
\end{equation} \tag{23} \]
где $n_p=11$ — число точек, по которым строится данная модель, то есть по которым находится оценка $\hat {c}_0$ аппроксимации временной зависимости факторной эластичности $p$ многочленом нулевой степени.

В третьем, пятом и последнем столбцах табл. 3 приведены оценки $s [p_k] = s [\eta _k ]$ среднеквадратичного отклонения результатов расчета $p_k$, $k=\overline {1,n_p}$. Так как очевидно, что эти результаты не являются равноточными, при оценке единственного параметра $c_0$ модели (23) можно воспользоваться взвешенным методом наименьших квадратов [12–15]. В качестве весов $\omega_k$ примем величины $\omega_k =s_k^2 = (s [p_k ] )^2$.

С учетом соотношения $\eta _k =\sqrt {\omega _k } \cdot \varepsilon _k$, где $\varepsilon _k $ — случайная величина, удовлетворяющая условиям теоремы Гаусса–Маркова [14], регрессию (23) преобразуем к виду
\[ \begin{equation}
\frac{p_k}{s_k}=\cfrac{1}{s_k}c_0 + \varepsilon_k, \quad k=\overline{1,np}.
\end{equation} \tag{24} \]
Заметим, что при этом из соотношения $s^2 [\eta _k ]=\omega _k s_\varepsilon ^2 =s_k^2$, где $s^2 [\eta _k ]$ и $s_\varepsilon^2$ — оценки дисперсии случайных величин $\eta _k$ и $\varepsilon _k$, следует, что $s_\varepsilon^2 = {s^2 [\eta _k ]/s_k^2=1}$, то есть при выбранных весах $\omega _k =s_k^2$ происходит нормирование оценки дисперсии случайного возмущения в регрессии (24).

Отсюда на основе минимизации $\| \varepsilon \|^2\to \min$ получаем оценку
\[ \begin{equation}
\hat {c}_0 = {\sum\limits_{k=1}^{n_p} p_k s_k^{-2}} \Big / {\sum\limits_{k=1}^{n_p} s_k^{-2}}.
\end{equation} \tag{25} \]
Оценку $s_{res}^2$ дисперсии $\sigma _\varepsilon^2 $ для модели (24) можно найти по формуле
\[ \begin{equation}
s_{res}^2 =\frac1{n_p-1} \sum _{k=1}^{n_p} s_k^{-2} ( p_k -\hat {c}_0 )^2.
\end{equation} \tag{26} \]

Во втором, третьем и четвертом столбцах табл. 4 приведены результаты расчетов оценок $\hat {c}_0$, $\sum _{k=1}^{n_p} s_k^{-2}$ и $s_{res}^2$ по описанным выше формулам (25) и (26) для факторных эластичностей $\alpha$, $\beta$ и $\gamma$ соответственно.

Таблица 4. Анализ зависимости факторных эластичностей от времени [Analysis of the dependence of factor elasticities on time]

Factor elasticities, $p$	$\hat {c}_0 \in \{ \hat {\alpha }_0,\hat {\beta }_0 ,\hat {\gamma }_0 \}$	$\sum\limits_{k=1}^{n_p} s_k^{-2}$	$s_{res}^2$	$s [c_0 ]$	$t$
$\alpha$	$-0.046$	1215.0	2.15	0.029	1.61
$\beta$	$-0.120$	1623.0	0.14	0.025	4.85
$\gamma$	$-0.645$	349.6	0.29	0.053	12.06

Проверку гипотезы об адекватности модели (24), то есть о корректности аппроксимации результатов расчета $p_k$, $k=\overline {1,n_p}$, описывающих временную зависимость факторной эластичности, константой $\hat {c}_0$, можно выполнить на основе распределения Фишера $F= s_{res}^2 / s_\varepsilon ^2$ [12]. При нормированной оценке дисперсии $s_\varepsilon ^2 =1$ имеем $F=s_{res}^2$ (четвертый столбец табл. 4). Сравнивая эти значения с $F_{cr} =F (0.05, 10, 18)=2.48$, имеем $F<F_{cr}$ для каждого параметра $\alpha$, $\beta$ и $\gamma$. Отсюда можно сделать вывод о правомерности выбора математической модели, описывающей динамику выпуска продукции региональной энергосистемы, в форме нелинейной трехфакторной степенной зависимости (1), факторные эластичности которой не меняются за период времени с 1990 по 2021 год.

Среднеквадратичное отклонение оценок $\hat {c}_0$, вычисленных по формуле (25) с учетом известных $s^2\big[p_k \big]=s_k^2$, можно найти на основе соотношений
\[ \begin{equation*}
s^2 [ \hat {c}_0 ]=\frac{ \sum _{k=1}^{n_p} ( s_k^{-2} )^2 s^2 [p_k ] }{\left( \sum _{k=1}^{n_p}
s_k^{-2} \right)^2}=
\frac{\sum _{k=1}^{n_p} s_k^{-2}}{\left(\sum _{k=1}^{n_p} s_k^{-2} \right)^2}
=
\frac{1}{\sum _{k=1}^{n_p} {s_k^{-2} } }.
\end{equation*} \]
Отсюда имеем, что
\[ \begin{equation}
s [\hat {c}_0 ]= \biggl( \sum\limits_{k=1}^{n_p} {s_k^{-2} } \biggr)^{-1/2}.
\end{equation} \tag{27} \]

Среднеквадратичные оценки $s [\hat {\alpha }_0 ]$, $s [ \hat {\beta }_0 ]$ и $s [ \hat {\gamma }_0 ]$, найденные по формуле (27), представлены в пятом столбце табл. 4. На их основе можно выдвинуть гипотезу о значимости параметра $\hat {c}_0 $ [12]. Для проверки этой гипотезы вычисляется статистика $t= | \hat {c}_0 |/ s [\hat {c}_0 ]$, значения которой для различных факторных эластичностей приведены в последнем столбце табл. 4. Сравнивая эти значения с критическим значением из таблицы распределения Стьюдента $t_{cr} =t (0.05, 10)=2.23$, можно сделать вывод о незначимости параметра $\alpha$ в нелинейной трехфакторной степенной зависимости (1). Это может служить поводом для упрощения модели (1), исключив из нее фактор $K( t )$, описывающий капитальные ресурсы. Однако из-за корреляции оценок параметров проверка значимости может оказаться не слишком точной, а упрощение модели — ненадежным. Как следствие, оценки прогноза получаются смещенными, и лучше оставить в модели незначимые параметры, чем отбросить значимые [12].

Таким образом, с учетом результатов проведенных численно-аналитических исследований и сравнительного анализа различных моделей функционирования региональной энергосистемы к практическому применению при решении задач оптимизации и планирования деятельности региональных энергопроизводств на 
основе системного анализа эффективности региональной энергетики можно рекомендовать ковариационно-стационарную модель временного ряда в форме 
разностного уравнения с мультипликативно-степенным трендом
\[ \begin{equation*}
y_k =\lambda _2 x_{1k}^{\lambda _3 } x_{2k}^{\lambda _4 } x_{3k}^{\lambda _5
} +\lambda _1 ( y_{k-1} -\lambda _2 x_{1,k-1}^{\lambda _3 } x_{2,k-1}^{\lambda _4 } x_{3,k-1}^{\lambda _5 } ) + \varepsilon _k
,\; k=1, 2, 3, \ldots, N-1.
\end{equation*} \]

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

Заключение

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

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

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

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

Конкурирующие интересы. Конкурирующих интересов не имеется.
Авторский вклад и ответственность. Все авторы принимали участие в разработке концепции статьи и в написании рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами.
Финансирование. Исследование выполнялось без финансирования.

About the authors

Vladimir E. Zoteev

Samara State Technical University

Email: zoteev.ve@samgtu.ru
ORCID iD: 0000-0001-7114-4894
SPIN-code: 8547-1223
Scopus Author ID: 16456013300
ResearcherId: D-8245-2014
http://www.mathnet.ru/person38585

Dr. Techn. Sci., Professor; Professor; Dept. of of Applied Mathematics and Computer Science

Russian Federation, 443100, Samara, Molodogvardeyskaya st., 244

Lyaysan A. Sagitova

Samara State Technical University

Email: l0410@mail.ru
ORCID iD: 0000-0002-0833-983X
SPIN-code: 5588-4106
https://www.mathnet.ru/person213377

Senior Lecturer; Dept. of Heat and Gas Supply and Ventilation

Russian Federation, 443100, Samara, Molodogvardeyskaya st., 244

Anna A Gavrilova

Samara State Technical University

Author for correspondence.
Email: a.a.gavrilova@mail.ru
ORCID iD: 0000-0001-6598-6518
https://www.mathnet.ru/person41413

Cand. Techn. Sci., Associate Professor; Associate Professor; Dept. of Control and System Analysis of Thermal Power and Sociotechnical Complexes

Russian Federation, 443100, Samara, Molodogvardeyskaya st., 244

References

Supplementary files