Стохастическое моделирование территориальных катастроф и преодоления их последствий Текст научной статьи по специальности «Математика»

СТОХАСТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ТЕРРИТОРИАЛЬНЫХ КАТАСТРОФ И ПРЕОДОЛЕНИЯ ИХ ПОСЛЕДСТВИЙ

Е.Н.Грибова

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

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

Ранее автором был составлен ряд программ, моделирующих кумулятивный стохастический рост с диссипацией, и высказаны идеи о различных типах управляющих воздействий, возможность внесения которых допускают эти программы [4, 6]. Некоторые из высказанных идей реализованы в данной работе.

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

m=input('ввод m:')

%m — масштаб квадрата S=rand(m);

%S — матрица случайных чисел m*m n=input('ввод n:')

%n — число шагов

^=трШ;('ввод времени (шага) взрыва tv:')

%tv — шаг наступления катастрофы k=1; p=0; vzr=0; g=0;

%счетчики while k<n for i=1:m

D(1,i)=rand; D(m,i)=rand;

D(i,1)=rand; D(i,m)=rand; end

for i=2:m-1 for j=2:m-1

D(i,j)=rand/4*(S(i-1,j)+S(i,j-1)+S(i+1,j)+S(i,j+1));

end

end

%нормируем приращение D NOR=sum(sum(D));

D(:,:)=0.5*mA2*D(:,:)/NOR;

%прирост уровня S(:,:)=S(:,:)+D(:,:); p=p+1; if p>9 for i=1:m for j=1:m D1(i,j)=rand;

%диссипация

S(i,j)=S(i,j)-2*D1(i,j);

end

end

p=0;

%рисунки

[X,Y]=meshgrid([1:1:m]); figure(1); plot3(X,Y,S) figure(2) meshc(X,Y,S) figure(3) surfc(X,Y,S) shading interp k

end

%цикл взрыва if k==tv

гу=трШ;('ввод координаты взрыва гу:') jv=input('ввод координаты взрыва _|у:') dv=input('ввод длины территории взрыва ^:') %if S(iv:(iv+dv),jv:(jv+dv))>FSr % if k==tv | vzr==1 %взрыв!!!

S(iv:(iv+dv),jv:(jv+dv))=rand;

'шаг взрыва' к

vzr=0; end if k>tv F=S; P=S;

F(iv:(iv+dv),jv:(jv+dv))=0;

%рамка

F(:,1:2)=0; F(:,(m-1):m)=0;

F(1:2,:)=0; F((m-1):m,:)=0;

'средние значения';

B=sum(sum(F));

FSr=B/(mA2-dvA2-8*m+16)

%F — поле без территории взрыва Р(:,1:^-1))=[]; P(:,(dv+2):(m-jv+1))=[]; P(1:(iv-1),:)=[]; Р^^^-^^М];

%Р — поле взрыва PSr=mean(mean(P))

RAT=PSr/FSr 'наибольшие значения';

Fmax=max(max(F));

Pmax=max(max(P));

к

if ЯАТ>0.5 'шаг восстановления' к

%печатаем шаг восстановления

vzr=1;

end

end

g=g+1; if g>9 R=S;

for z=1:12

%z — количество расчетов max во всей матрице [rs,ir]=max(R);

Wmax=max(rs);

for i=1:m

if Wmax==rs(i)

a(z)=ir(i); b(z)=i;

end

end

STR=struct('nom',z,Wmax',Wmax,'i',a(z),'j',b(z))

R(a(z),b(z))=0;

end

g=0;

ttt=input('пауза')

end

k=k+1;

end

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

Периодически (а именно, в данной программе — через каждые 10 шагов цикла) происходит случайное «осыпание» полученного рельефа, при котором значение уровня каждой клетки уменьшается на учетверенную (при «щадящем» режиме диссипации — удвоенную) случайную величину, лежащую в тех же пределах.

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

Периодически (в представленном варианте программы — через каждые 10 шагов цикла) на печать выводятся рисунки, иллюстрирующие текущее состояние рельефа, а также данные о некотором количестве (в данной версии — о двенадцати) наиболее высоких точек этого рельефа: уровень этих точек и их координаты в системе 0, j). Тем самым, данная программа позволяет проследить формирование центра (элемента рассматриваемой территориальной системы с максимальной высотой) и бицентра (второго по уровню элемента системы) [5], стохастический характер их становления, развития, а в некоторых случаях и упадка.

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

На заранее задаваемом шаге ^ программа моделирует локальную катастрофу (условно говоря — взрыв) на территории квадрата с вершиной (гу, jv) и длиной стороны dv, значения которых вводятся с клавиатуры непосредственно на ^-м шаге исполнения программы. «Взрыв» заключается в том, что каждой точке выбранного квадрата приписывается значение, равное случайной величине,

лежащей между 0 и 1. Эта же функция, естественно, позволяет осуществить и «перезапуск» всего процесса («глобальный» взрыв), если на ^-м шаге выбирается квадрат с вершиной (1, 1) и длиной стороны т.

Программа предусматривает возможность осуществления многократных «взрывов» на различных участках исходного поля. Кроме того, имеется возможность программировать взрыв в зависимости от достижения определенных условий — например, взрывать всегда только бицентр либо предусматривать эпицентр взрыва в точке с определенным заранее заданным уровнем, и т.д.

Рис. 1. Моделирование катастрофы центра в условиях кумулятивного роста (т=50, п=10)

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

величин RAT=PSr/FSr. В случае наступления определенных условий, задаваемых исследователем (определенное соотношение максимальных по модулю или средних элементов поля взрыва и остальной части территории), программа выдает сообщение о факте восстановления поврежденной взрывом территории.

Рис. 2. Моделирование катастрофы центра в условиях кумулятивного роста (т=50, п=40)

Характер работы данной программы при наступлении катастрофы показывают рисунки 1 и 2, иллюстрирующие последовательные состояния одного и того же процесса. На 11-м шаге произошел взрыв квадрата 4х4 с вершиной в точке (13, 23). «Сдвоенный» центр территории, который отчетливо виден на рисунке 1, таким образом, оказался разрушен, и его место занял бывший бицентр, к 40-му шагу разросшийся в целую группу точек. Помимо него, на рисунке 2 хорошо различима группа точек, составляющая бицентр, и изолированная вершина 9 с координатами (8, 29).

Рисунки 3 и 4 иллюстрируют ситуацию, когда ни центр, ни бицентр непосредственно не пострадали от взрыва (tv = 18, iv = 27, jv = 23, dv = 6), однако, в силу кумулятивного характера процесса роста, группа точек, сформировавшая на 10-м шаге бицентр (28,

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

Картина с ликвидацией последствий катастрофы выглядит достаточно благополучно (к 40-му шагу RAT неизменно достигает 0,35-0,45), если программа допускает кумулятивный стохастический рост без насыщения. Ситуация претерпевает коренные изменения в случае, когда совокупный объем распределяемых на каждом шаге ресурсов ограничен определенной величиной, возможно, зависящей от времени (это некий аналог распределенного ВНП). Задавая различную зависимость ВНП от времени, можно в случае необходимости моделировать динамику экономического подъема, кризиса, депрессии, циклические колебания конъюнктуры на данной территории.

Рис. 3. Моделирование становления бнцентров в условиях катастрофы (т=50, п=20)

Рис. 4. Моделирование становления бицентров в условиях катастрофы (т=50, п=40)

Рис. 5. Моделирование катастрофы в условиях кумулятивного роста с насыщением (т=50, п=30)

m=50, п=70

С

£

а эс

Рис. 6. Моделирование катастрофы в условиях кумулятивного роста с насыщением (т=50, п=70)

Катастрофу в условиях кумулятивного роста с насыщением и диссипацией иллюстрируют рисунки 5 и 6, отражающие соответственно 30-й и 70-й шаги одного и того же процесса. Несмотря на то, что первоначально поле взрыва начало в среднем подниматься (Р8г выросло с -0,6713 на 20-м шаге до -0,5404 на 30-м), затем началось стабильное ухудшение: Р8г(40) = -3,1678, Р8г(50) = -6,6712, Р8г(60) = -12,9881, Р8г(70) = -17,9789. Обращает на себя внимание тот факт, что средняя глубина территории взрыва увеличилась с 20го по 70-й шаг более чем в 25 раз, тогда как высота центра, не изменившего за этот период своих координат, выросла лишь в 4 раза — с 20,6046 на 20-м шаге до 81,5786 на 70-м.

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

«щадящей» диссипации (т=50, п=30)

Рис. 8. Моделирование катастрофы в режиме

«щадящей» диссипации (т=50, п=60)

Даже введение «щадящего» режима диссипации, проиллюстрированного рисунками 7 и 8, отражающими различные стадии одного и того же процесса, позволяет лишь смягчить, но не устранить до конца указанную общую закономерность (PSr(30) = -

0,4902, PSr(50) = -1,9322, PSr(70) = -4,1629). При этом следует отметить тот факт, что территория, потерпевшая катастрофу, нередко имеет тенденцию к расширению, захватывая и соседние участки поля, на которых тоже начинается обвальное падение уровня развития, как это видно на рисунках 7 и 8 [7].

Радикальное преодоление отмеченной тенденции возможно в случае, когда часть территории взрыва (даже одна точка) начинает развиваться в режиме «ударной стройки», т.е. получает «целевым назначением» некоторый объем ресурсов за счет пропорционального замедления роста остальных точек всего поля. Управляющее воздействие такого рода порождает различные сценарии быстрого выхода из кризиса части территории, потерпевшей катастрофу, что хорошо согласуется с многочисленными реальными примерами регионального развития, в частности, нашей страны (освоение целинных земель, восстановление Ташкента после землетрясения 1966 года и проч.).

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

ш=три1;('ввод m:')

%т-масштаб квадрата S=rand(m);

%S-матрица случайных чисел m*m п=трШ;('ввод n:')

%п-число шагов %счетчики:

%шаг; шаг взрыва; рост; восстановление; k=1; kv=100; р=0; vzr=0;

%статистика; взрыв; стройка; g=0; tv=0; ss=0; ts=0;

for i=1:m

D(1,i)=rand; D(m,i)=rand;

D(i,1)=rand; D(i,m)=rand; end while k<n if k>2

D(1,1)=rand/4*(S(1,2)+S(2,1)); D(m,m)=rand/4*(S(m,m-1)+S(m-1,m)); D(1,m)=rand/4*(S(1,m-1)+S(2,m)); D(m,1)=rand/4*(S(m-1,1)+S(m,2)); for i=2:m-1

D(1,i)=rand/6*(S(1,i-1)+S(2,i)+S(1,i+1)); D(m,i)=rand/6*(S(m,i-1)+S(m-1,i)+S(m,i+1)); D(i,1)=rand/6*(S(i-1,1)+S(i,2)+S(i+1,1)); D(i,m)=rand/6*(S(i-1,m)+S(i,m- 1)+S(i+1,m)); end end

if k<kv | k==kv | (k>kv & RAT>1) | ss==0 for i=2:m-1 for j=2:m-1

D(i,j)=rand/4*(S(i-1,j)+S(i,j-1)+S(i+1,j)+S(i,j+1));

end

end

%нормируем приращение D NOR=sum(sum(D));

D(:,:)=0.5*mA2*D(:,:)/N0R;

%прирост клеток

S(:,:)=S(:,:)+D(:,:);

P=P+1;

end

%стройка

if k>kv & RAT<1 & ss==1 %по всему for i=2:m-1 for j=2:m-1

D(i,j)=rand/4*(S(i-1,j)+S(i,j-1)+S(i+1,j)+S(i,j+1));

end

end

%по стройке %по запросу поля for i=(is-ds):(is+ds) for j =(j s-ds):(j s+ds) D(i,j)=rand*(S(i-1,j)+S(i,j-1)+S(i+1,j)+S(i,j+1)); end end

N0R=sum(sum(D));

D(:,:)=0.5*mA2*D(:,:)/N0R;

%прирост клеток S(:,:)=S(:,:)+D(:,:);

P=P+1;

end

%осыпание if p>9 for i=1:m for j=1:m D1(i,j)=rand;

S(i,j)=S(i,j)-4*D1(i,j);

end

end

P=0;

%рисунки

[X,Y]=meshgrid([1:1:m]);

figure(1);

Plot3(X,Y,S) figure(2) surfc(X,Y,S) shading interP figure(3) surfc(X,Y,S) k end

%цикл взрыва if tv==1

1у=три1;('ввод координаты взрыва гу:') jv=input('ввод координаты взрыва _|у:') dv=input('ввод длины территории взрыва dv:')

%if S(iv:(iv+dv),jv:(jv+dv))>FSr k==tv | vzr==1 %взрыв!!!

S(iv:(iv+dv),jv:(j v+dv))=rand;

'шаг взрыва' kv=k

vzr=0; tv=0; end if ^==1

88=трЩ;('строшъ сейчас?^) да=1, нет=0:') if Б8==1

is=input('ввод координаты стройки is:') js=input('ввод координаты стройки js:') ds=input('ввод длины территории стройки ds:') end ts=0; end

if k==kv | к>к^

F=S; P=S;

F(iv:(iv+dv),jv:(jv+dv))=0;

%рамка

F(:,1:2)=0; F(:,(m-1):m)=0;

F(1:2,:)=0; F((m-1):m,:)=0;

'среднее значение';

B=sum(sum(F));

FSr=B/(mA2-dvA2-8*m+16)

%F-поле без территории взрыва Р(:,1:^-1))=[]; P(:,(dv+2):(m-jv+1))=[]; P(1:(iv-1),:)=[]; Р^^^-^^М]; %Р-поле взрыва PSr=mean(mean(P))

RAT=PSr/FSr 'наибольшее значение';

Fmax=max(max(F));

Pmax=max(max(P)); if RAT>0.5 'шаг восст' k

%печатаем шаг восстановления vzr=1; end end

g=g+1;

%статистика

if g>9 R=S;

for z=1:12

%z-количество расчетов мах во всей матрице [rs,ir]=max(R);

Wmax=max(rs); for i=1:m if Wmax==rs(i) a(z)=ir(i); b(z)=i; end end

STR=struct('nom',z,Wmax',Wmax,'i',a(z),'j',b(z))

R(a(z),b(z))=0;

end

g=0;

tv=input('взрывать?(tv) да=1:') ts=input('строить?(ts) да=1:') end k=k+1; end

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

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

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

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

Рис. 9. Моделирование катастроф и преодоления их последствий в условиях кумулятивного роста с насыщением (т=50, п=10)

Рис. 10. Моделирование катастроф и преодоления их последствий в условиях кумулятивного роста с насыщением (т=50, п=20)

Рис. 11. Моделирование катастроф и преодоления их последствий в условиях кумулятивного роста с насыщением (т=50, п=30)

Однако, если на этой же территории происходит повторный взрыв (21-й шаг процесса), то восстановление проходит значительно дольше (вплоть до 45-го шага) и характеризуется нелинейной стохастической динамикой (рис. 11).

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

Рисунки 12-15 являются иллюстрациями различных состояний одного и того же процесса.

Рис. 12. Моделирование катастроф и преодоления их последствий в условиях кумулятивного роста с насыщением (т=50, п=10)

На рис. 12 отчетливо видны «сдвоенные» центр и бицентр.

На 11-м шаге произошел взрыв квадрата 4х4 с вершиной в точке (26, 37) — центр территории, который оказался разрушен.

Рис. 13. Моделирование катастроф и преодоления их последствий в условиях кумулятивного роста с насыщением

(т=50, п=20)

При пассивном наблюдении (второй этап динамики рассматриваемого региона) территория центра медленно, но постепенно начинает восстанавливаться. Но после «планового» осыпания (шаг 20) его уровень значительно снизился (рис. 13). Кроме того, определились новый, опять “сдвоенный”, центр (24, 46) и бицентр (39,

42).

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

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

Рис. 14. Моделирование катастроф и преодоления их последствий в условиях кумулятивного роста с насыщением (т=50, п=30)

На четвертом этапе при введении стройки по всей территории взрыва рост параметров идет еще более быстрыми темпами и к 47му шагу происходит восстановление поврежденных территорий. Стабилизируется депрессивный характер развития донорских районов — как непосредственно прилежащих к разрушенной территории, так и бывшего бицентра, послужившего в данном приме-

ре основным донором. При этом новым бицентром данного региона становится соседний с прежним бицентром район.

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

Рис. 15. Моделирование катастроф и преодоления их последствий в условиях кумулятивного роста с насыщением (т=50, п=50)

Несмотря на небольшое осыпание нового бицентра и попытку введения дополнительной стройки на территории районов-доноров, в короткие сроки изменить ситуацию не удается (рис. 15). Стоит заметить, что быстро растущий новый бицентр к 50-му шагу стабильно закрепил за собой статус центра, а бицентром соответственно становится территория, выступавшая до этого центром.

Представленный комплекс имитационных стохастических моделей позволяет: 1) отслеживать неравномерность роста различных составных частей замкнутых территорий, 2) фиксировать закономерности территориальной динамики, присутствующие в развитии отдельных стран и регионов, 3) моделировать в различных режи-

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

Литература

1. Нижегородцев Р.М. Траектории региональной экономики: проблемы причинности и моделирования // Информация и экономика: теория, модели, технологии / Под ред. Е.Ю.Иванова, Р.М.Нижегородцева. Барнаул: АГУ, 2002. С. 172-194.

2. Нижегородцев Р.М. Информационная экономика. Книга 3. Взгляд в Зазеркалье: Технико-экономическая динамика кризисной экономики России. Москва-Кострома, 2002. — 170 с.

3. Нижегородцев Р.М. Поляризация экономического пространства России и как ей противодействовать//Проблемы теории и практики управления. 2003. № 1. — С. 89-95.

4. Нижегородцев Р.М., Грибова Е.Н. Апология бицентризма: рост и упадок, анализ и прогнозирование// Региональная экономика в информационном измерении: модели, оценки, прогнозы: Сб. науч. трудов/Под ред. Е.Ю.Иванова, Р.М. Нижегородцева. Москва — Барнаул: Изд-во «Бизнес-Юнитек», 2003. — С. 81-114.

5. Нижегородцев Р.М., Грибова Е.Н. Моделирование возникнове-

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

практической конференции. Новомосковск, 2003. — С. 47-63.

6. Грибова Е.Н., Нижегородцев Р.М. Управление развитием территорий: имитационное моделирование катастроф и преодоления их последствий//Управление большими системами. Вып. 5. М.: ИПУ РАН, 2003. — С. 21-34.

7. Грибова Е.Н., Нижегородцев Р.М. Моделирование локальных траекторий регионального развития: катастрофы и преодоление их последствий// Проблемы управления безопасностью сложных систем. Часть 1. М.: РГГУ, 2003. С. 147-151.