Автоматика и телемеханика, № 8, 2019

Управление в социально-экономических

системах

(Государственный университет штата Миннесота, Сен-Клауд, США)

ДВОЙСТВЕННЫЙ АЛГОРИТМ ПРОГНОЗИРОВАНИЯ

ТЕХНОЛОГИЧЕСКИХ МАТРИЧНЫХ СТРУКТУР

В ДИНАМИЧЕСКИХ МОДЕЛЯХ ТИПА ЗАТРАТЫ—ВЫПУСК

На основе глобального метода последовательного улучшения Крото-

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

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

че квадратичного программирования большой размерности с сепарабель-

ным функционалом, возникающей при прогнозировании матрицы пря-

мых затрат в динамических моделях межотраслевого баланса. С помо-

щью декомпозиции удается воспользоваться специальным видом матри-

цы ограничений для уменьшения размерности задачи.

Ключевые слова: модель межотраслевого баланса, матрица прямых за-

трат, сбалансированное прогнозирование, квадратичное программиро-

вание, декомпозиция, двойственный метод оптимального управления

В.Ф. Кротова.

DOI: 10.1134/S0005231019080099

1. Введение

Для эффективного принятия решений и прогнозирования на макроэко-

номическом уровне необходимы информационные и программные средства,

обеспечивающие анализ данных, имитационное моделирование и расчет оп-

тимальных сценариев на основе математических моделей и методов системно-

го анализа. В современной экономике сфера применения таких компьютер-

ных систем включает разработку макроэкономических сценариев развития

на уровне страны (региона) с учетом межотраслевого баланса.

Основы моделей межотраслевого баланса (моделей типа затраты—выпуск)

были заложены В.В. Леонтьевым [1] и развиты многочисленными последо-

вателями. В работах В.И. Гурмана, В.Ф. Кротова и их коллективов [2-6]

были предложены обобщения для модели Леонтьева, в том числе были по-

строены оптимизационные модели для социо-эколого-экономических систем.

Привлечение аппарата теории оптимального управления динамическими си-

стемами [7-9] стало основой для создания технологий для расчета сценариев

в динамике с учетом тех или иных управлений, целевых критериев качества

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

сматриваемые в таких моделях, нередко оказываются достаточно сложны-

ми в математическом и вычислительном смыслах (например, ввиду возмож-

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

109

Л.С. Понтрягина, большой размерности системы и других факторов), то для

решения таких задач оказывается не достаточным применение тех или иных

«универсальных» методов теории оптимального управления и соответствую-

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

числительных алгоритмов и программных средств, учитывающих специфику

моделей.

В частности, в работах коллективов В.И. Гурмана и В.Ф. Кротова (в Ин-

ституте проблем управления РАН, и Институте динамики систем и теории

управления Сибирского отделения РАН) было создано специальное матема-

тическое и программное обеспечение для идентификации систем и построе-

ния динамических сценариев на основе управлений (например, [2-6, 10, 11]).

Соответствующие методы решения задач оптимального управления во мно-

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

ных В.Ф. Кротовым [9].

Цель статьи — разработка вычислительного алгоритма, основанного на

известном методе глобального улучшения В.Ф. Кротова, для сбалансиро-

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

затраты—выпуск, описанной ниже, с учетом ее специфики.

2. Постановка задачи

Рассматривается динамическая модель следующего общего вида [12]:

∑

(2.1)

Xⁱ(t) =

a_ij(t)X^j(t) + Yⁱ(t) - Impⁱ

(t),

j=1

∑

(2.2)

Yⁱ(t) =

d_ij(t)I^j(t) + Cⁱ(t) + Expⁱ

(t),

j=1

(2.3)

Kⁱ(t + 1) = Kⁱ(t) + Vⁱ(t) - Wⁱ

(t),

(2.4)

V i(t) = ΦiV (Iⁱ(t),... ,Iⁱ(t - τvi

)),

(2.5)

Wⁱ(t) = Φ^iW (Kⁱ(t),... ,Kⁱ(t - τwi

)),

(2.6)

Xⁱ(t) = Fⁱ(t,Kⁱ(t),Lⁱ

(t)),

(2.7)

A(t) = Φ_A(X(t), . . . , {A(tkA

)}),

(2.8)

D(t) = Φ_D(I(t), . . . , {D(tkD

)}),

где i = 1, n, t = t₀, T (например, годы), Xⁱ(t) — валовые выпуски; Yⁱ(t) — ко-

нечный продукт; Impⁱ(t), Expⁱ(t) — соответственно импорт и экспорт i-й от-

расли; Cⁱ(t) — непроизводственное потребление продукции отрасли i (вклю-

чающее в себя личное потребление, расходы на оборону, прочие статьи го-

сударственного потребления, а также запасы); Kⁱ(t) — стоимость основных

производственных фондов (ОПФ); Vⁱ(t), Wⁱ(t) — величина ввода в действие

и выбытие ОПФ отрасли i; Iⁱ(t) — капитальные вложения в i-ю отрасль в

году t; Lⁱ(t) — численность занятых в i-й отрасли; A(t) = ∥a_ij(t)∥ — матрица

110

коэффициентов прямых затрат; {A(tkA )}, {D(tkA )} — множества отчетных

статистических матриц; D(t) = ∥d_ij (t)∥ — матрица технологической структу-

ры капитальных вложений. Модель включает в себя балансовые соотношения

блоков с постоянной структурой (2.1)-(2.3) и блоков с переменной структу-

рой (2.4)-(2.8), для каждого из которых может существовать более чем один

вариант соответствующих им операторов F_i, Φ^iV , Φ^iW , Φ_A, Φ_D, i = 1, n. От-

метим, что многие практически интересные производственные функции при-

дают модели нелинейный характер.

В общем случае данные соотношения представлены в виде уравнений с

распределенными лагами (2.4) и (2.5), где τvi и τwi обозначают максималь-

ную глубину запаздывания соответственно в процессах обновления и старе-

ния фондов i-й отрасли.

В данной статье акцент сделан на один из центральных блоков модели —

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

ства, т. е. прогнозирование матричной функции A(t).

С внедрением в производство новых технологий происходит переход на

использование новых материалов, изменение их пропорций, что влечет изме-

нение количественных показателей, характеризующих межотраслевые связи.

Применительно к модели это означает, что изменению подвергается матрица

коэффициентов прямых затрат.

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

способами. Один из них — задать тенденции изменения коэффициентов пря-

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

технологических изменений, т. е. выяснить, как такие изменения повлияют

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

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

связывающие технологические коэффициенты с остальными переменными

модели.

Ряд методов прогнозирования динамики матрицы A(t) предполагает ста-

тистическое экстраполирование тенденций в изменении ее элементов. Такие

алгоритмы (см., например, [13, 14]) позволяют достаточно точно отразить

поведение коэффициента каждого в отдельности, но не гарантируют адек-

ватного описания общих (суммарных) тенденций в материальных затратах и

промежуточном потреблении.

Достижение сбалансированности матричных элементов с их суммами по

строкам и столбцам (строчными и столбцовыми окаймлениями) обеспечива-

ет известный метод RAS [15]. Основным предположением в нем является то,

что более устойчивые и гладкие траектории, присущие окаймлениям, долж-

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

метод полностью абстрагируется от изменений в матрице под влиянием вре-

менного фактора и в явном виде время вообще не учитывает.

Модернизированные методы прогнозирования матриц A(t), D(t) синте-

зируют в себе оба подхода и предназначены для получения «структурно-

динамического сбалансированного прогноза» [16]. Идея такого подхода со-

стоит в нахождении прогнозных значений элементов матрицы, динамика ко-

торых определяется поведением их траекторий в предшествующий период и

111

X_ĳ

P_i

Z_j

Рис. 1. Окаймления матрицы межотраслевых потоков.

факторами, воздействующими на них в рассматриваемый период. В отличие

от метода RAS в данном методе фактор времени учитывается в явном виде.

Рассматриваемый ниже метод, основанный на [16], применительно к прогно-

зированию матрицы прямых затрат был развит в [12].

Изложим суть метода в обобщенном виде применительно к матрице пря-

мых затрат A. Введем обозначения: X_ij(t) — матрица межотраслевых пото-

∑

ков; p_i =

X_ij — вектор промежуточного продукта отраслей, идущего на

j=1

∑

производственное потребление; z_j = X_ij — вектор материальных затрат

i=1

∑

отраслей; Q = p_i = z_j =

∑X_ij — совокупная величина промежуточ-

i j

ного продукта. Таким образом, p_i и z_j , являющиеся соответственно строч-

ными и столбцовыми окаймлениями матрицы X_ij , сбалансированы, в свою

очередь, с величиной Q (рис. 1).

Аналогично для матрицы коэффициентов прямых затрат a_ij = X_ij /X^j

имеет место сбалансированность:

∑

(2.9)

p_i =

a_ijX^j, z_j =

a_ijX^j

j = 1,n, i = 1,m.

j=1

i=1

Если ввести величины α_i = p_i/Xⁱ — долю промежуточного продукта в ва-

ловом выпуске отраслей, β_j = z_j /X^j — материалоемкость валового продукта

∑

отраслей, x_i = Xⁱ/ X^j — удельный вес i-й отрасли в валовом внутреннем

j=1

продукте (ВВП), то тождество (2.9) можно переписать в виде:

∑

(2.10)

β_j(t) =

a_ij

(t),

i=1

∑

(2.11)

α_i(t) =

a_ij(t)x_j(t)/x_i

(t).

j=1

112

∑ n∑

Величина γ =

X_ij/ X^j

= Q/X, интерпретируемая: а) как доля

i j

суммарного промежуточного продукта в ВВП и б) как материалоемкость

общественного продукта, соответственно равна

∑

(2.12)

γ(t) =

α_i(t)x_i(t) =

β_j(t)x_j

(t),

i=1

j=1

где β_j (t), α_i(t), γ(t) — соответственно элементы окаймляющей строки, окайм-

ляющего столбца и суммарная величина окаймления, x_j(t), x_i(t) — в общем

случае весовые множители, суть коэффициенты отраслевой структуры вало-

вого продукта.

В предлагаемом методе выделяются два подхода — дезагрегационный

и комбинированный [12], отличающиеся способом согласования прогнозных

значений окаймляющих величин с прогнозными величинами â_ij(t).

При дезагрегационном подходе алгоритм прогнозирования разбивается на

три этапа.

1. Осуществляется прогноз суммарной доли γ: ее расчетные величины γ(t)

формируются либо с помощью подбора кривых роста, либо экспертным пу-

тем.

2. Строятся “частные” (не зависимые друг от друга) прогнозы α_i(t)

β_j(t)

окаймляющих долей α_i и β_j, а затем производится их согласование (баланси-

ровка) с суммарной долей γ(t) в соответствии с (2.12). Для этого решаются

две задачи квадратичного программирования [17]

⎧

∑

⎪

⁽α_i(t) - α_i(t))2

⎪

-→ min,

⎪

α_i(t)

⎨

i=1

(2.13)

∑

⎪

α_i(t)x_i(t) = γ(t),

⎪

⎩ i⁼¹

α_i(t) ≤ α_i(t) ≤ α_i(t), i = 1,m,

⎧

(

)₂

⎪

∑

β_j(t)

⎪

-→ min,

⎪

⎨

β_j(t)

j=1

(2.14)

∑

⎪

β_j(t)x_j(t) = γ(t),

⎪

⎩ j⁼¹

β_j(t) ≤ β_j(t)

β_j(t), j = 1,n,

в которых минимизируется суммарное относительное отклонение искомых

величин α_i(t) и β_j (t) от их частных прогнозов. При этом учитываются ниж-

ние и верхние границы допустимого варьирования переменных α_i(t) и β_j (t),

соответственно α_i(t), α_i(t) и β_j(t),

β_j(t). Отметим, что прогноз весовых ко-

эффициентов x_i(t) также осуществляется автономно.

Таким образом, при решении задач (2.13), (2.14) согласуются три группы

независимо сформированных прогнозов: x_i(t), γ(t) и {α_i(t)

β_j(t)} (рис. 2).

113

a_ĳ

X ⁱ

X ^j

Рис. 2. Согласование независимых прогнозов.

3. Прогноз коэффициентов прямых затрат â_ij(t) осуществляется на основе

модификации предложенного в [16] “динамизированного” метода RAS.

Предполагается, что поведение расчетных траекторий элементов данной

матрицы â_ij(t) определяется близостью их к имеющимся в отдельные годы

отчетного периода элементам матриц A(t₁), A(t₂), . . ., A(t_k) и изменением

элементов базисной матрицы A⁰ под воздействием прогнозных окаймлений.

(Не нарушая общности рассуждений, можно за базисную матрицу A⁰ по-

ложить последнюю, имеющуюся в статистической отчетности A(t_k), тогда

начальный год моделирования t₀ = t_k). При этом предполагается, что “отно-

сительный прирост любого элемента структурной матрицы в году t форми-

руется из двух аддитивных составляющих, одна из которых является общей

для всех элементов данного столбца, другая — общей для всех элементов

данной строки” [16]. Такая зависимость отражается формулой

(

)

(2.15)

â_ij(t) = a0ij

1 + r_if_i(t) + s_jg_j(t)

Функции f_i(t) и g_j (t), такие что f_i(t₀) = g_j (t₀) = 0, i = 1, m, j = 1, n, опре-

деляют характер изменения составляющих относительных темпов прироста

во времени. Их вид считается известным, и характер поведения в основном

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

Прогноз коэффициентов матрицы â_ij(t) с учетом сделанных выше предпо-

ложений о характере поведения расчетных траекторий сводится к необходи-

мости решения задачи квадратичного программирования (ЗКП) следующего

вида:

⎧

⎪

∑[

(

)]₂

⎪

a0ij

1 + r_if_i(t_l) + s_jg_j(t_l) - a_ij(t_l)

-→ min,

⎪

r_i,s_j

⎪

i=1 j=1 l=1

⎪

∑

(

)

⎨

a0ij

1 + r_if_i(t) + s_jg_j(t)

x_j = α_i(t)x_i, i = 1,m,

(2.16)

⎪

j=1

⎪

∑

(

)

⎪

a0ij

1 + r_if_i(t) + s_jg_j(t)

= β_j(t)), j = 1,n,

⎪

⎩ i⁼¹

(

)

a_ij(t) ≤ a0ij

1 + r_if_i(t) + s_jg_j(t)

≤ a_ij(t), i = 1,m, j = 1,n.

114

Статистические

Последовательность

Ограничения,

и эндогенные

дезагрегации

экзогенно

переменные

показателей

задаваемые

промежуточного

пользователем

потребления

x_i(t)

(t)

i(t)

j(t)

i(t), j(t)

i(t)

j(t)

f_i(t), g_j(t)

r_i, s_j

ĳ(t)

{

(t)

ĳ(t )},

ĳ(t)

l = 1, k

Рис. 3. Схема дезагрегационного алгоритма.

В данной постановке в отличие от [16], где рассматривается прогноз пото-

ков, а не коэффициентов, имеют место также дополнительные ограничения

типа неравенств, отражающие точку зрения экономистов, строящих прогноз.

С помощью величин a_ij (t), a_ij(t) можно управлять поведением отдельных,

наиболее важных коэффициентов a_ij(t).

На схеме (рис. 3) изобразим причинно-следственные связи в описанном

дезагрегационном алгоритме прогнозирования.

Характерной особенностью “дезагрегационного” прогнозирования являет-

ся то, что в нем неявно предполагают, что достоверность прогноза суммар-

ной величины выше, чем достоверность суммы прогнозов ее составляющих

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

вероятнее всего имеют в виду свойство наколения вычислительной ошибки

при суммировании и то, что динамика окаймлений имеет более устойчивый

и гладкий вид по сравнению с динамикой суммированных величин. Однако

именно последний факт позволяет рассчитывать на то, что в данном случае

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

имно компенсироваться. Поэтому согласование расчетных траекторий â_ij(t) с

прогнозами α_i(t), β_j (t) (предварительно согласованными с γ(t)) имеет смысл

не только для повышения надежности прогноза матрицы A(t), сколько по-

115

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

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

Таким образом, в процессе дезагрегирования информации по цепочке

(γ) → (α_i, β_j ) → (a_ij) матрица A(t) “теряет” статус параметра модели и ста-

новится эндогенной переменной, причем не промежуточной, а итоговой. Дей-

ствительно, во-первых, из выражений (2.14)-(2.16) следует зависимость A(t)

от x_i, а главное, во-вторых, если матрица A(t) рассчитана в соответствии с

системой (2.16) при уже известных величинах α_i(t), то для нахождения век-

тора конечного продукта {Yⁱ} нет необходимости решать уравнения меж-

отраслевого баланса Y = (E - A) X, поскольку тогда достаточно вычислить

Yⁱ = (1 - α_i)Xⁱ, i = 1,m. Другими словами, при прогнозировании коэффи-

циентов прямых затрат в рамках данного подхода речь идет не об “учете

динамики параметра A(t)”, а о прогнозе собственно межотраслевых связей,

моделируемых технологическими коэффициентами a_ij .

Алгоритм прогнозирования при комбинированном подходе заключается в

независимом построении прогнозов γ(t), α_i(t)

β_j(t) и их одновременном со-

гласовании с теоретической зависимостью (2.15) для a_ij(t, r_i, s_j). При этом

в отличие от описанной выше трехэтапной дезагрегационной процедуры ре-

шается одна, комбинированная задача, в которой минимизируется средне-

взвешенное квадратическое отклонение â_ij (t_l) от их статистических вели-

чин a_ij (t_l) и искомых γ, α_i, β_j от их частных прогнозов γ(t), α_i(t)

β_j(t):

∑

∑(

)₂

θ_a

â_ij(t_l,r_i,s_j) - a_ij(t_l)

i=1 j=1 l=1

(

)₂

∑

⁽α_i - α_i(t))2

β_j

β_j(t)

(2.17)

+θ_α

+θ_β

α_i(t)

β_j(t)

i=1

j=1

⁽γ - γ(t))2

+θ_γ

-→ min

γ(t)

r_i,s_j,α_i,β_j,γ

∑

(2.18)

â_ij(t,r_i,s_j)x_j(t) = α_ix_i

(t), i = 1, m,

j=1

∑

(2.19)

â_ij(t,r_i,s_j) = β_j

j = 1,n,

i=1

∑

(2.20)

â_ij(t,r_i,s_j)x_j

(t) = γ,

i=1 j=1

(2.21)

a_ij(t) ≤ â_ij(t,r_i,s_j) ≤ a_ij

(t),

i = 1,m, j = 1,n,

(2.22)

α_i(t) ≤ α_i(t) ≤ α_i

(t),

(2.23)

β_j(t) ≤ β_j(t)

β_j

(t),

(2.24)

γ(t) ≤ γ(t) ≤ γ(t).

116

Здесь θ_a, θ_α, θ_β, θγ — задаваемые экзогенно весовые множители, с по-

мощью которых регулируется суммарная степень близости каж(дой из)груп-

пы показателей к соответствующим статистическим значениям

a_ij(t_l)

либо

(

)

частным прогнозам

γ(t), α_i(t)

β_j(t)

В данной задаче ограничение (2.20) совместно с (2.18), (2.19) обеспечивает

выполнение балансовых равенств (2.12) и поэтому может быть эквивалентно

заменено на любое из них.

Отметим, что в “комбинированном” алгоритме прогноз матрицы A(t)

влияет на расчет конечного продукта Yⁱ = (1 - α_i)Xⁱ, поскольку в зависи-

мости (2.17)-(2.24) величины γ, α_i, β_j , â_ij оцениваются совместно в отличие

от алгоритма (2.14)-(2.16).

Численные методы прогнозирования для описанных алгоритмов предла-

гаются в следующем разделе.

Следует отметить, что в рамках природно-экономической модели [4] этот

метод можно применить также для прогнозирования матрицы коэффициен-

тов прямых затрат на восстановление ресурсов A(z) и матрицы коэффициен-

тов фондообразующих затрат при восстановлении ресурсов B(z).

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

технологических параметров модели

Рассмотренные выше алгоритмы прогнозирования межотраслевых мат-

ричных структур, дезагрегационный и комбинированный, предполагают ре-

шение для каждого года моделирования t = [t₀, T ] задач квадратичного про-

граммирования (2.13), (2.14), (2.16) либо (2.17)-(2.24).

В силу большой размерности решение задач (2.16) и (2.17)-(2.24) извест-

ными методами нелинейного программирования представляется весьма за-

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

может оказаться эффективным применение двойственного итерационного ме-

тода [7]. Данный алгоритм является универсальным методом решения задачи

отыскания закона управления и соответствующих ему фазовых траекторий

таких, чтобы процесс не выходил бы за пределы заданных ограничений, с

одновременной минимизацией некоторого критерия. Он опирается на доста-

точные условия оптимальности управляемых процессов, которые являются

эффективными в том случае, если может быть найдена некоторая функция ϕ

состояния процесса, обладающая рядом специальных свойств.

3.1. Двойственный метод оптимизации многошаговых процессов

Прежде чем представить вычислительный алгоритм, напомним основные

конструкции двойственного метода глобального улучшения В.Ф. Кротова

[7, 9]. Для сохранения терминологии статьи [7] в данном пункте принятые

ранее обозначения временно переопределяются.

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

ментом t, состоянием y и управлением u, причем t = 0, N - 1, y ∈ Y = Rⁿ,

117

u ∈ U ⊂ R^m, и удовлетворяющий уравнению

(3.1)

y(t + 1) = f(t, y(t), u(t)), y(0) = y₀, y_i(N) = b_i

i = 1,r ≤ n.

Такие процессы будем называть допустимыми, а их множество обозначим

через D. Определим на множестве E = {(y(t), u(t)) : y^Rn, u ∈ U} функционал

вида

∑

(3.2)

J = f⁰

(t, y(t), u(t)) + F (N, y).

t=0

Требуется найти последовательность {y_s(t), u_s(t)} элементов из E, сходящую-

ся к D и минимизирующую J на D:

(3.3)

J (y_s, u_s) → d = inf

Определим в соответствии с [7] класс функций ϕ(t, y) ∈ Π и следующие

конструкции:

(3.4)

R(t, y, u) = ϕ(t + 1, f(t, y, u)) - ϕ(t, u) - f⁰

(t, y, u),

(3.5)

G(N, y) = ϕ(N, y) + F (N, y),

∑

(3.6)

L = G(N,y) - R(t,y,u) - ϕ(0,y₀

t=0

∑

(3.7)

l = min G(N,y) -

maxR(t, x, u) - ϕ(0, y0

yr+1,...,y

y,u

t=0

Справедливы следующие утверждения [7].

Лемма 1. ∀ ϕ ∈ Π:

1) L(y, u) = J(y, u), если (y, u) ∈ D,

2) l(ϕ) ≤ d.

Теорема 1. Если (y,u) ∈ D и ϕ ∈ Π такие, что L(ϕ,y,u) = l(ϕ), то

J (y, u) = min J(y, u) ≡ d = l(ϕ) = max l(ϕ).

Процесс (y, u), удовлетворяющий условиям теоремы, называется опти-

мальным. Для его отыскания будем строить последовательность {ϕ_s} ∈ Π

начиная с производной ϕ₀, на которой l(ϕ_s) → max l(ϕ).

Будем искать ϕs+1 в виде

(3.8)

ϕs+1(t, y) = ϕs(t, y) + λsγs(t, y),

где λ_s > 0, γ_s(t, y) - подлежащие определению скаляр и функция.

Определим функционал вида

∑

[

]

(3.9)

δ_s(y,u) = γ_s(N,b) -

γ_s(t + 1,f(t,y,u)) - γ_s(t,y)

- γ_s(0,y₀

t=0

118

который можно переписать следующим образом:

(3.10)

∑[

]

∑



δ_s(y,u) =

γ_s(t,y) - γ_s(t,y(t) - z(t - 1))

γ_s(t,y)

t=tk⁺⁰

t=t_k -0

t=1

t₁=0, t₂=N

где z(t) = y(t + 1) - f(t, y, u).

Для разности Δl = l_s(λ) - l_s можно показать, что

(3.11)

Δl_s(λ) = λδ_s(y(t),u(t)) + {L_s(y(t),u(t)) - l_s

где (y(t), u(t)) ∈ E_s(λ). Из уравнений (3.6), (3.7) следует неотрицательность

второго слагаемого в формуле (3.11).

Элементарная операция улучшения функции ϕ_s(t, y) состоит в построении

ϕs+1(t, y) в виде зависимости (3.8), где λs > 0 и γs выбираются произвольно,

так чтобы

(3.12)

δ_s

(y, u) > 0

хотя бы при одном значении (y, u) ∈ Es+1 = E_s(λ_s).

Из выражения (3.11) непосредственно следует, что элементарная опера-

ция (3.8) удовлетворяет неравенству

(3.13)

Δl_s = ls+1 - l_s

> 0,

поэтому может использоваться для построения улучшающей последователь-

ности {(y_s, u - s)} : l_s → d.

3.2. Применение двойственного глобального метода

к решению выпуклой задачи квадратичного программирования

с сепарабельным функционалом

С учетом специфики данных задач для решения их описанным выше ите-

рационным методом [7] необходимо провести их декомпозицию и представить

в виде дискретной задачи оптимального управления вида (3.1), (3.2).

Рассмотрим ЗКП следующего вида [12]:

∑

(3.14)

J (u) =

(c1hu2h + c0hu_h) -→ min,

h=1

∑

(3.15)

c_ihu_h = b_i

i = 1,n,

h=1

(3.16)

uminh ≤ u_h ≤ umaxh

h = 1,N,

∥C_ih∥, b_i, C1h > 0, C0h, uminh, umaxh — заданные величины.

119

Она является выпуклой задачей квадратичного программирования с сепа-

рабельным функционалом¹. В данном случае функционал представлен мно-

гочленом.

Задачу (3.14)-(3.16) можно представить в форме эквивалентной задачи

оптимизации многошагового процесса с одномерным аргументом t = 0, N,

скалярным управлением u_t и состоянием x_t ∈ Rⁿ:

∑

(3.17)

J (u) =

(c1tu2t + c0tu_t) -→ min,

t=1

(3.18)

x_i(t + 1) = x_i(t) + c_i(t)u_t

i = 1,n, t = 0,N - 1,

(3.19)

x_i

(0) = 0, i = 1, n,

x_i(N) = b_i, i = 1,n,

{

}

(3.20)

u_t ∈ U =

u_t : umint ≤ u_t ≤ umaxt

Для решения поставленной задачи воспользуемся двойственным методом

оптимизации многошаговых процессов.

Зададим функцию u(t, x) в линейном виде

∑

(3.21)

ϕ(t, x) =

ψ_i(t)x_i

i=1

и определим для нее конструкции (3.24) и (3.25)

∑[

]

R(t, x, u) =

ψ_i(t + 1) - ψ_i(t)

x_i+

i=1

(3.22)

[

]

∑

+ ψ_i(t + 1)c_i(t) - c⁰

u_t - c1tu2t,

i=1

∑

(3.23)

G(N, x) = ψ_i(N) x_i = ψ_i(N) b_i.

i=1

Пояснение:

(

)

R(t, x, u) = ψ

t + 1,f(t,x,u)

- u(t, x) - f⁰(t, x, u) =

∑

[

]

∑

[

]

= ψ_i(t + 1)

x_i + c_i(t)u(t)

- ψ_i(t)x_i -

c1tu2t + c0tu_t

i=1

[

]

∑

= xi [ψ_i(t + 1) - ψ_i(t)] +

ψ_i(t + 1)c_i(t) - c⁰

- u²(t)c1t.

i=1

¹ Многочлен f(x) называется сепарабельным, если его неприводимые множители не

имеют кратных корней.

120

Условие существования максимума по x и функции R(t, x, u) накладывает

на ϕ ∈ Π следующие условия (3.24), (3.26):

(3.24)

ψ_i(t + 1) - ψ_i(t) = 0, т.е. ψ_i(t + 1) = ψ_i(t) - ψ_i

= const.

Пояснение:

Данное положение вытекает из необходимого условия максимума

∂R

= ψ_i(t + 1) - ψ_i(t) = 0

∀i = 1,n,

∂x_i

∑

∂R

(3.25)

= -2c1tu + ψ_i(t + 1)c_i(t) - c0t

= 0.

∂u

i=1

Таким образом функция R не зависит от x:

⎧

⎨

umint,

ũ≤umint,

u_t =

ũ_t,

umint < ũ_t < umaxt,

⎩

umaxt,

ũ_t ≥ umaxt,

(3.26)

ⁿ ψ_ic_i(t) - c0t

ũ_t =i=1

2c¹

Это также вытекает из необходимого условия максимума (см. выше).

Двойственный функционал (3.27) будет следующего вида:

∑

{[

]

}

(3.27)

l(ψ) =

ψ_ib_i -

ψ_ic_i(t) - c0t

u_t(ψ) - c1tu2t(ψ)

i=1

t=0

Ставится задача:

⎧

⎪

l(ψ) -→ max,

⎨

∑

⎪

⎩ c

ihuh = bi, i = 1, n.

h=1

Функцию γ будем задавать в соответствии с (3.21) и (3.23):

∑

(3.28)

γ(t, x) =

ν_ix_i.

i=1

Тогда функция δ_s примет вид

(

)

∑

(3.29)

δ_s(x,u) =

ν_i b_i -

c_i(t) ũ_t

i=1

t=0

121

Пояснение:

∑[

(

)]

∑

δ_s(x,u) =

γ_s(t,x) - γ_s

t,x(t) - z(t - 1)

+ γ_s(t,x) |t=tk+0t=t

k-0

t=1

t₁=0

∑

[

(

)

]

= γ_s(N,b) -

γ_s

t + 1,f(t,x,u)

- γ_s(t,x)

- γ_s(0,x₀) =

t=0{

(

)

(

)

= γ_s(N,b) - γ_s

1, f(0, x, u)

- γ(0,x) + γ_s

2, f(1, x, u)

}

(

)

-γ_s(1,x) + ... + γ_s

N,f(N - 1,x,u)

- γ_s(N - 1,x)

- γ_s(0,x₀)-

−{ ... } = +γ_s(0,x) - γ_s(N,x) =

= γ_s(N,b) - γ_s(N,x) + γ_s(0,x) - γ_s(0,x₀) =

∑

= ν_ib_i - ν_i c_i(t) u_t + 0 - ν_i0 =

i=1

t=0

i=1

(

)

∑

= ν_i b_i - c_i(t) u_t

i=1

t=0

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

обеспечения δ_s > 0, положим

∑

(3.30)

ν_i = b_i -

c_i(t) u_t.

t=0

Таким образом, можем записать следующий алгоритм построения мини-

мизирующей последовательности:

1 шаг: Задаемся произвольными значениями ψ0i .

2 шаг: ψs+1i = ψ^si + λ_sνsi , где

λ_s = argmax ls+1(λ)

λ>0

ls+1(λ) = l(ψs+1) = l(ψ^s + λν^s).

3 шаг: Если ls+1 - l_s < ε, то выход из алгоритма, иначе переход ко второму

шагу.

3.3. Построение численного алгоритма прогнозирования.

Декомпозиция задачи

Декомпозиционными или блочными принято называть такие методы ре-

шения оптимизационной задачи, которые связаны с анализом составляющих

ее подзадач [17].

Членение на подзадачи осуществляется до начала процесса решения, при-

чем, как правило, выбор метода декомпозиции существенно зависит от этого

122

членения. Метод членения обычно реализуется в виде последовательности

итераций, на каждой из них производится независимый анализ отдельных

подзадач при фиксированных значениях параметров, связывающих подзада-

чи. От итерации к итерации значения этих параметров пересчитываются [18].

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

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

дача охватывает все ограничения и часть переменных общей задачи, во вто-

ром — имеет место двойственная ситуация: отдельные блоки включают все

переменные и некоторые ее ограничения. Также существуют методы, которые

допускают одновременное использование вертикального и горизонтального

членения задачи, как например, метод, предложенный в [19].

В целевых функциях алгоритмов прогнозирования (2.16) и (2.17)-(2.24)

матрицы квадратичной формы помимо диагональных членов содержат так-

же коэффициенты при смешанных произведениях r_i · s_j. Для того, чтобы из-

бавиться от них и привести задачу к ЗКП с сепарабельным функционалом,

проведем декомпозицию: будем попеременн{ фиксировать оди} из векторов

r_i и s_j и решать последовательность задач

ЗКП1p(r), ЗКП2p(s)

1. Шаг 0. Задаем произвольно u0j, j = 1, N ;

ψ⁰}.

2. Шаг р. Решаем последовательно две задачи

⎧

⎪ЗКП1p : управление — upi = ri, sj = up-1j = fixe, ls -→max, ψp0

ψp-1;

⎨

ψ^s

⎪

⎩ЗКП^p

: управление — up+1j = s_j, r_j = u^pi = fixe, ls+1 -→

max, ψp+10

ψ^p.

ψs+1

3. l

ψp+1) - l

ψp-1) < ϵ, то выход; иначе переход к п. 2.

Несмотря на то, что максимум двойственного функционала l^p в каждой

из задач ЗКП1p, ЗКП2p не обеспечивает строгого выполнения ограничений —

равенств (уравнений процесса), для сходимости описанной и{ерационной про-

цедуры достаточно, чтобы для любого p выполнялось U^p =

u^p : u^min(up-1) ≤

}

≤u^p ≤u^max(up-1)

= 0. Тогда будет иметь место улучшение l^p(ψ), ∀p, что в

силу построения алгоритма и обеспечивает сходимость.

Далее рассмотрим применение данной итеративной декомпозиционной схе-

мы для решения задач в данной модели (очевидно, что задачи (2.13) и (2.14)

тривиально записываются в виде (3.14)-(3.16) и не требуют декомпозиции):

Задача (2.16)

ЗКП1p: Запишем задачу в терминах постановки (3.14)-(3.16):

u_h = r_i (h = i = 1,N), s_j - fixe (N = n - число отраслей).

∑∑[

]₂

(3.31)

c1h =

a0hjf_h(t_l)

h = 1,N;

j=1 l=1

∑∑[

]

(3.32)

c0h = 2

a0hj(1 + s_jg_j(t_l)) - a_hj(t_l)

a0hjf_j(t_l

h = 1,N;

j=1 l=1

123

⎧⎛

⎞

⎪

∑

⎨

⎝ a⁰

⎠fh(t)x_h · δ_ih, i = 1,N (δ_ih - символ Кронекера);

(3.33)

c_ih =

j=1

⎪

⎩a0hi-N · f_h(t), i = N + 1,2N, h = 1,N;

⎧

∑

⎪

α_ix_i -

a0ij(1 + s_jg_j(t))x_j, i = 1,N;

⎨

(3.34)

b_i =

[

]

⎪

∑

(

)

⎪

⎩ βi-N -

a0m

1 + s_i-Ng_i-N(t)

i = N + 1,2N;

i-N

{

}

ahj(t) - ahj · (1 + sjgj(t))

(3.35)

uminh = max

, h = 1,N;

j,f_h(t)=0

f_h(t)

{

}

ahj(t) - ahj · (1 + sjgj(t))

(3.36)

umaxh = min

, h = 1,N.

j,f_h(t)=0

f_h(t)

ЗКП2p: u_h = s_j (h = j = 1, N), r_i - fixe.

∑∑[

]₂

(3.37)

c1h =

a0ihg_h(t_l)

h = 1,N;

i=1 l=1

∑∑[

]

(3.38)

c0h = 2

a0ih(1 + r_if_i(t_l)) - a_ih(t_l)

a0ihg_h(t_l

h = 1,N;

i=1 l=1

⎧

⎨a0ihgh(t)xh, i = 1,N;

(

)

(3.39)

c_ih =

∑

⎩

a⁰

g_h(t) · δ_i-N,h, i = N + 1,2N, h = 1,N;

m=1

⎧

⎡

⎤

⎪

∑

⎪

α_ix_i -⎣ a⁰

x_j⎦(1 + r_if_i(t)), i = 1,N;

⎨

j=1

(3.40)

b_i =

⎪

∑

(

)

⎪

⎩β_i-N - a0m,i-N

1 + r_mf_m(t)

i = N + 1,2N;

m=1

{a_ih(t) - a0ih · (1 + r_if_i(t))},h=1,N;

(3.41)

uminh = max

i,g_h(t)=0

g_h(t)

}

{a_ih(t) - a0ih · (1 + r_if_i(t))

(3.42)

umaxh = min

, h = 1,N.

i,g_h(t)=0

g_h(t)

Задача (2.17)-(2.24)

Декомпозиция данной задачи также основана на попеременном фиксиро-

вании одного из векторов r_i, s_j. Однако управление u_h здесь имеет для каж-

124

ЗКП

i = 1, n

i = n + 1, 2n

i = 2n + 1

Рис. 4. Представление задач в соответствии с компонентами управления u_h =

= (r, α, β, γ), h = 1, N.

дой частной задачи не n, a N = 3n + 1 компонент: ЗКП1p : u_h = {r_i, v_i, w_j , u},

ЗКП2p : u_h = {s_j, v_i, w_j , u}.

Структура матрицы системы линейных ограничений c_ih имеет вид, пока-

занный на рис. 4, где изображен вариант, когда вместо ограничения (2.22)

∑_n

записано эквивалентное уравнение

v_ix_i = u (n — общее число отраслей

i=1

в модели).

Формулы для коэффициентов c1h, c0h, c_ih, b_i, α_h, β_h в данной задаче в отно-

шении первой составляющей переменной управления u_h: r_i или s_j совпадают

с приведенными выше выражениями.

Остальные элементы перечисленных векторов и матрицы коэффициентов,

соответствующие компонентам v_i, w_j , u управления u_h одинаковы для обеих

задач ЗКП1p, ЗКП2p и могут быть легко записаны непосредственно исходя из

формулировки задачи (2.17)-(2.24).

3.4. Реализация и функционирование алгоритма

Описанный выше декомпозиционный алгоритм инвариантен относительно

вида теоретической зависимости для â_ij(t) в рамках класса линейных двух-

параметрических функций времени вида

(3.43)

â_ij = σ1ijf_ij(t) + σ0ijg_ij

(t), i, j = 1, n.

Таким образом, помимо рассмотренной аддитивной формулы (2.15) мож-

но использовать различные кривые роста типа (3.43), применяемые обыч-

но в одномерном прогнозировании. Несмотря на существенное увеличение

при этом количества оцениваемых параметров (вместо 2n переменных r_i,

s_j - 2n² переменных σ1ij, σ0ij), можно предположить, что быстродействие ал-

горитма практически не ухудшится. Это следует из того, что при декомпози-

ции σ1ij либо σ0ij становятся значениями управляющей переменной u_h при h =

= i(n - 1) + j; h = 1, N; N = n² + 2n + 1, а сложность задачи (3.17)-(3.18)

зависит не от N, а от размерности фазового вектора y, т.е. от общего числа

ограничений-равенств. Если также учесть в реализации алгоритма сильную

разреженность матрицы c_ih, то можно сократить вычислительные затраты

при суммировании по h = 1, N в несколько раз.

Что касается способов конкретизации зависимости (2.15), то функции f_i(t)

и g_j(t), которые должны отражать динамику окаймлений [16], можно непо-

125

средственно задавать в виде относительных темпов прироста соответствую-

щих долей v_i, и w_j :

v_i(t₁) - v_i(t₀)

w_j(t₁) - w_j(t₀)

(3.44)

f_i(t₁) =

, g_j(t₁) =

, i, j = 1, n, l = 1, k.

v_i(t₀)

w_j(t₀)

В задаче (2.17)-(2.24) для t₁ = t величины v_i, w_j в (3.44) следует заменить

на их частные прогнозы v_i(t), ŵj(t), û(t) и априорных ограничениий vi(t),

vi(t), wj (t), wj(t), u(t), u(t). В системе предусмотрены два способа.

1. Частные прогнозы определяются из рекуррентных формул

(

)

(

)

vi(t) =

1+δvi

v_i(t - 1),

ŵ_j(t) =

1+δwj

w_j(t - 1), û(t) = (1 + δ_u)u(t - 1),

где темпы прироста δvi , δwj , δ_u задаются пользователем, так же как и дина-

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

несколько видов кривых роста.

2. Частные прогнозы строятся с помощью подбора временных зависимо-

стей из заданного набора путем оценивания параметров по рядам статисти-

ческих наблюдений. Найденные параметры затем пользователь может кор-

ректировать для задания требуемых тенденций в прогнозном периоде. Огра-

ничения в данном случае формируются в виде «трубки»:

(

)

(

)

1-Δvi

≤v_i ≤v_i

1+Δvi

(

)

(

)

ŵ_j

1-Δwj

≤w_j ≤ ŵ_j

1+Δwj

û(1 - Δ_u) ≤ u ≤ û(1 + Δ_u)

с заданными величинами допустимых относительных отклонений Δvi , Δwj ,

иΔ_u.

Аналогично можно двумя способами формировать ограничения на a_ij.

1. Строятся частные прогнозыâij(t), например, в виде уравнения (3.43), и

тогда

a_ij(t) = (1 - Δ_ij)âij(t),

a_ij(t) = (1 + Δ_ij)âij(t), i,j = 1,n.

При решении задач согласования происходит уточнение параметров σ1ij,

σ0ij зависимости (3.43). Отметим, что для эффективного построения функ-

ций (3.43) требуется представительная динамика отчетных матриц {a_ij (t₁)},

l = 1,K, где K ≥ 15 - 20.

2. В случае отсутствия длинных статистических рядов {a_ij (t)} можно

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

(

)

â_ij(t) =

1+δaij

a_ij(t - 1) и строить «трубку» ограничений с помощью до-

пустимых отклонений Δ_ij.

Для наиболее важных коэффициентов прямых затрат (их доля не превы-

шает 20%) можно жестко фиксировать плановые значения a_ij(t), т. е. Δ_ij = 0.

В процессе функционирования алгоритма пользователь помимо указанных

экзогенных величин может управлять параметрами работы алгоритма θ_a, θ_y,

θ_w, θ_u, {ε_j} для достижения требуемой точности и обеспечения сходимости и

совместности решаемой задачи.

126

4. Заключение

На основе метода глобального улучшения В.Ф. Кротова в работе пред-

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

цы A(t) прямых затрат в рассмотренной динамической модели межотрасле-

вого баланса.

В описанных процедурах прогнозирования сочетается возможность

экспертно-аналитического формирования частных прогнозов показателей

промежуточного потребления с учетом требования их сбалансированности

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

позволяет моделировать достаточно широкий спектр технологических изме-

нений в структуре материального производства.

Декомпозиция алгоритмов прогнозирования и применение для решения

задач квадратичного программирования (ЗКП) с выпуклым сепарабельным

функционалом двойственного метода оптимизации многошаговых процес-

сов [7], основанного на операции улучшения функции Кротова, дает возмож-

ность существенно увеличить быстродействие алгоритмов по сравнению с ис-

пользованием стандартных программ решения ЗКП общего вида.

Исследования, изложенные в данной статье, были начаты во время работы

автора в Институте проблем управления (ИПУ) РАН имени В.А. Трапезни-

кова в лаборатории 45 под руководством В.Ф. Кротова. Автор признателен

сотрудникам лаборатории А.Г. Александрову, О.В. Моржину и Л.А. Селива-

новой за обсуждение результатов и помощь в работе.

СПИСОК ЛИТЕРАТУРЫ

1. Leontief W.W. Input-Output Economics. 2nd Ed. N.Y.: Oxford Univer. Press, 1986.

2. Оптимальное управление природно-экономическими системами / Под ред.

В.И. Гурмана, А.И. Москаленко. Новосибирск: Наука, 1980.

3. Кротов В.Ф. Исследование нелинейных оптимизационных моделей развития

многоотраслевой экономики. Ч. I-III // АиТ. 1981. № 10. C. 129-136; № 11.

C. 114-123; 1982. № 1. C. 114-122.

Krotov V.F. Investigation of Nonlinear Optimization Models of Growth in a

Multisectoral Economy. I-III // Autom. Remote Control. 1981. V. 42. No. 10.

P. 1385-1391; No. 11. P. 1524-1531; 1982. V. 43. No. 1. P. 91-98.

4. Эколого-экономическая стратегия развития региона / Под ред. В.И. Гурмана.

Новосибирск: Наука, 1990.

5. Моделирование социо-эколого-экономической системы региона / Под ред.

В.И. Гурмана, Е.В. Рюминой. М.: Наука, 2003.

6. Proops J., Safonov P. Modeling in Ecological Economics / Edward Elgar Publ.

U.K., 2004.

7. Кротов В.Ф. Вычислительные алгоритмы решения и оптимизации управляе-

мых систем уравнений I, II // Изв. АН СССР. Техн. кибернетика. 1975. № 5,

№ 6.

8. Кротов В.Ф., Фельдман Н.Н. Итерационный метод решения задач оптималь-

ного управления // Изв. АН СССР. Техн. кибернетика. 1983. № 2. С. 160—168.

9. Krotov V.F. Global Methods in Optimal Control Theory. N.Y.: Marcel Dekker,

1996.

127

10. Гурман В.И., Матвеев Г.А., Трушкова Е.А. Социо-эколого-экономическая мо-

дель региона в параллельных вычислениях / Управление большими системами.

Сб. тр. Вып. 32. М.: ИПУ РАН, 2011. С. 109-130.

11. Расина И.В., Блинов А.О., Гусева И.С. Магистрали в задаче оптимизации стра-

тегии развития региона на многокомпонентной модели // Вест. Бурят. гос. ун-

та. Сер. 9. Математика и информатика. 2011. С. 36-42.

12. Сафонов П.И. Алгоритм прогнозирования прямых затрат в динамической

модели межотраслевого баланса / Вопросы создания АСПР. Применение

экономико-математических методов в перспективном планировании. Сб. научн.

трудов ГВЦ Госплана СССР. М.: 1988. Вып. 87. С. 118-138.

13. Николаева И.Г., Новикова Т.А. Метод статистического прогнозирования пря-

мых затрат / Прикладные задачи экономического моделирования. Вып. 13. М.:

ВНИИСИ, 1984.

14. Николаева И.Г. Анализ и прогноз межотраслевых связей. М.: Экономика, 1981.

15. Stone R., Brown J.A.C. A Long-Term Growth Model of the British Economy //

Eur. Future Figures. Amsterdam, 1962.

16. Седелев Б.В., Журавлева Л.В. Сбалансированный прогноз структурно-дина-

мических компонент макроэкономических процессов. Прикладные задачи эко-

номического моделирования / Сб. тр. ВНИИСИ. М., 1984. Вып. 13. С. 87-95.

17. Ляшенко И.Н. и др. Агрегирование и декомпозиция в моделях народнохозяй-

ственного планирования. Киев: Знание, 1980.

18. Овсеенко О.А., Шемякина Т.Ю. Методы декомпозиции процессов управления.

М.: МИУ, 1986.

19. Гольштейн Е.Г. Метод декомпозиции задач линейного и выпуклого програм-

мирования / Экономика и мат. методы. М., 1985. Вып. 6. Т. 21.

Статья представлена к публикации членом редколлегии М.М. Хрусталевым.

Поступила в редакцию 14.02.2017

После доработки 25.07.2018

Принята к публикации 08.11.2018

128